“Code for:
Ecology and sexual conflict drive the macroevolutionary dynamics of female-limited colour polymorphisms”
Beatriz Willink
2025-06-12
This documents describes the analyses conducted by Willink et al. (2024) to investigate the macroevolutionary dynamics of female-limited colour polymorphisms (FP) in pond damselflies (Coenagrionoidae).
Load packages
x <-
c(
"tidyverse",
"ggplot2",
"wesanderson",
"coda",
"ape",
"convenience",
"data.table",
"geiger",
"ggtree",
"gridExtra",
"treeio",
"phangorn",
"phytools",
"RevGadgets",
"kableExtra",
"MCMCglmm",
"MASS"
)
lapply(x, function(y) {
# check if installed, if not install
if (!y %in% installed.packages()[, "Package"])
install.packages(y)
# load package
try(require(y, character.only = T), silent = T)
})Data preparation
Tree pruning
We based our analysis on the MAP tree of Willink et al. (2024). This tree includes both pond damselflies (Coenagrionidae) and featherlegs (Platycnemididae). I’m not aware of any instance of FP in featherlegs after searching all species in Willink et al. (2024). So for this study we’ll focus only on pond damselflies.
I start by pruning the tree in Willink et al. (2024) to leave only pond damselfly taxa with sex-related colour data.
# read the tree
tre <- read.nexus("../data/phylogeny/G1_beta_strong.673464_MAP.tree")
# read character data
femdat <- read.csv("../data/processed/HiSSE/Fem_colour.csv", header = T, sep = ",")
# filter missing data
femdat <- femdat[femdat$State != "-",]
# check names in data and tree match
missing_taxa <- name.check(tre, femdat, femdat$Taxon)
# drop tips with missing character data
tre <- ape::drop.tip(tre, missing_taxa$tree_not_data)
# get featherleg ancestor
featherleg <- findMRCA(tre, c("Elattoneura_caesia", "Arabicnemis_caerulea"))
# drop featherlegs from tree
tre <- ape::drop.tip(tre, Descendants(tre, featherleg, type = "tips")[[1]])
# filter featherlegs from character data
feather_taxa <- name.check(tre, femdat, femdat$Taxon)
femdat <- femdat[!(femdat$Taxon %in% feather_taxa$data_not_tree),]
# write pruned tree
#write.nexus(tre, file = "../data/processed/HiSSE/Willink_MAP_coen.tree", translate = F)
# write filtered character data
#write.table(femdat, file = "../data/processed/HiSSE/Fem_colour_coen.csv", quote = F, row.names = F, sep = "\t")Summary of the character state data
For starters, I’d like to know how many species we have in each category.
# read the final data
sta <- read.table(file = "../data/processed/HiSSE/Fem_colour_coen.csv", header = T, sep = "\t")
sta_summ <-
data.frame(
state = c("SM", "SD", "FP", "missing", "Total"),
ntaxa = c(
length(sta$State[sta$State == "0"]),
length(sta$State[sta$State == "2"]),
length(sta$State[sta$State == "1"]),
length(sta$State[sta$State == "-"]),
length(sta$State)
)
)
sta_summ %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| state | ntaxa |
|---|---|
| SM | 116 |
| SD | 183 |
| FP | 119 |
| missing | 0 |
| Total | 418 |
Transition rates and diversification
We used a Hidden-State Speciation and Extinction (HiSSE) model to simultaneously infer transition rates between sex-related colour states and the effects of sexual dimorphism (SD) and female-limited colour polymorphisms (FP) on diversification dynamics.
Running the model
The model was implemented on RevBayes v 1.1 (Höhna et al. 2016). I ran two independent chains of this analysis. I did preliminary runs on UPPMAX. I ran the final model on my laptop, which has RevBayes installed via spack.
Processing the output
Convergence assessments
I used the convenience package to assess stationarity and convergence. While some ESS values fall under the 625 threshold in individual runs, the lowest from an individual run is 328 and the lowest ESS from the joint runs is 964. Moreover all parameter posteriors converge between the two chains.
check_HiSSE <- checkConvergence( list_files = c("../output/HiSSE/PolyHiSSE_MAP_coen_r1.model.log", "../output/HiSSE/PolyHiSSE_MAP_coen_r2.model.log"), control = makeControl(burnin=.10, precision = 0.01))## Reading in log file 1
## Reading in log file 1
## [1] "Analyzing continuous parameters"## FAILED CONVERGENCE
##
## 15 parameters failed to reach 625 for ESS_run_1
## 5 parameters failed to reach 625 for ESS_run_2
##
## BURN-IN SET AT 0.1
##
##
## LOWEST CONTINUOUS PARAMETER ESS
## RUN 1 -> speciation.3. 328.33
## RUN 2 -> speciation_observed.3. 606.35
## ## means ESS
## extinction.1. 0.0013911741 7744.2308
## extinction.2. 0.0468767168 3529.9953
## extinction.3. 0.0015391372 7983.0000
## extinction.4. 0.0036436526 7563.9857
## extinction.5. 0.0974009925 3239.1282
## extinction.6. 0.0072090397 7075.3084
## extinction_hidden.1. 0.7041097449 6635.8913
## extinction_hidden.2. 1.2958902593 6635.8921
## extinction_hidden_sd 0.5135588920 6480.5421
## extinction_hidden_unormalized.1. 0.7426905346 6632.1868
## extinction_hidden_unormalized.2. 1.5514405850 7339.7085
## extinction_observed.1. 0.0025174133 7621.4143
## extinction_observed.2. 0.0721388532 3098.5198
## extinction_observed.3. 0.0043740876 7050.5928
## hidden_rate 0.0002794787 1676.6248
## R.1. 0.0002794787 1676.6242
## R.2. 0.0002794787 1676.6242
## speciation.1. 0.0321656469 1805.3599
## speciation.2. 0.0971110810 3236.0433
## speciation.3. 0.0525591856 1142.2080
## speciation.4. 0.1416434600 1099.8813
## speciation.5. 0.4082545735 1192.5535
## speciation.6. 0.2317871038 1009.6954
## speciation_hidden.1. 0.4163896420 1049.0333
## speciation_hidden.2. 1.5836103658 1049.0320
## speciation_hidden_sd 1.0335103967 1287.6001
## speciation_hidden_unormalized.1. 0.5116576957 1058.5675
## speciation_hidden_unormalized.2. 2.0614595002 2208.9916
## speciation_observed.1. 0.0869045551 1083.0266
## speciation_observed.2. 0.2526828270 1284.8029
## speciation_observed.3. 0.1421731388 964.7533
Run under prior
I ran the same model without data, i.e. under the prior
#!bin/bash
#load RevBayes
source ~/Applications/spack/share/setup-env.sh
spack load revbayes /sum2kc5
# run diversification analysis on pruned Coenagrionidae tree
rb_command="source(\"./HiSSE/HiSSE_setup_MAP_prior_pruned.Rev\");"
echo $rb_command | rbAnd quickly confirm that posterior transition and diversification rates match the priors
# read posterior of run under prior
pst <- read.table("~/Dropbox/Doctorado/COEN_Comparative/FP_Comparative/output/HiSSE/PolyHiSSE_MAP_coen_prior.model.log", header = T)
# transitions
# reshape to long
prior_t <- pivot_longer(pst, cols = c(42:47), names_to = "Transition", values_to = "Rate")
# get colour palette
pal <- wes_palette("Zissou1", n = 6, type = "continuous")
# plot rates
p_t <- ggplot(prior_t,aes(x = Rate, fill = Transition)) +
theme_bw(base_size = 8) +
labs(x = "Transition rate", y = "Posterior frequency") +
geom_histogram(
binwidth = 0.0001,
linewidth = 0.5,
alpha = 0.6,
colour = "white",
position = "identity"
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
#xlim(0, 1.0) +
scale_fill_manual(
values = pal) +
theme_minimal(base_size = 7) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
# speciation
# reshape to long
prior_s <- pivot_longer(pst, cols = c(28:33), names_to = "Speciation", values_to = "Rate")
# plot rates
p_s <- ggplot(prior_s,aes(x = Rate, fill = Speciation)) +
theme_bw(base_size = 8) +
labs(x = "Speciation rate", y = "Posterior frequency") +
geom_histogram(
bins=100,
linewidth = 0.5,
alpha = 0.6,
colour = "white",
position = "identity"
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
xlim(-1, 99.0) +
scale_fill_manual(
values = pal) +
theme_minimal(base_size = 7) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
# extinction
# reshape to long
prior_e <- pivot_longer(pst, cols = c(5:10), names_to = "Extinction", values_to = "Rate")
# plot rates
p_e <- ggplot(prior_e,aes(x = Rate, fill = Extinction)) +
theme_bw(base_size = 8) +
labs(x = "Extinction rate", y = "Posterior frequency") +
geom_histogram(
bins=100,
linewidth = 0.5,
alpha = 0.6,
colour = "white",
position = "identity"
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
xlim(-1, 99.0) +
scale_fill_manual(
values = pal) +
theme_minimal(base_size = 7) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
lay_matrix <- rbind(c(1,1),
c(2,3))
grid.arrange(p_t, p_s, p_e, layout_matrix = lay_matrix)
Ancestral states
First, let’s look at the ancestral state estimation. This is to get a sense of the history of sex-related colouration in pond damselflies, but not to address any specific hypothesis.
Run the annotation on RevBayes
#!bin/bash
#load RevBayes
source ~/Applications/spack/share/setup-env.sh
spack load revbayes /sum2kc5
rb_command="source(\"./HiSSE/Annot_tree.Rev\");"
echo $rb_command | rb# read the annotated tree - it doesn't matter which run we use
stree_fn = "../output/HiSSE/PolyHiSSE_MAP_coen_r1_ase.tre"
stree_fn = "../output/HiSSE/PolyHiSSE_MAP_coen_r2_ase.tre"
# process ancestral states, here we focus on the colour trait only
tre <- processAncStates(stree_fn,
state_labels = c("0" = "SM_1",
"1" = "FP_1",
"2" = "SD_1",
"3" = "SM_2",
"4" = "FP_2",
"5" = "SD_2"))## | | | 0% | |========================================| 100%# define colour palette
st_colors = c( wes_palette("Zissou1", n = 12, type = "continuous")[c(12,12,7,7,1,1)], "white" )
# plot tree with ancestral states as pies
zz=plotAncStatesPie(t=tre,
pie_colors=st_colors,
node_labels_size=0,
node_pie_size = 0.75,
#node_pie_nudge=0.5,
tip_pie_size=0.5,
tip_pie_nudge_x = 1,
node_pie_nudge_x = 0.5,
tip_labels=FALSE,
xlim=c(0,106),
state_transparency = 0.75,
layout="circular"
)
zz
For completeness, we can plot the ancestral states for the hidden trait that influences diversification dynamics.
# process ancestral states, here we focus on the colour trait only
tre <- processAncStates(stree_fn,
state_labels = c("0" = "slow_FP",
"1" = "slow_SM",
"2" = "slow_SD",
"3" = "fast_FP",
"4" = "fast_SM",
"5" = "fast_SD"))## | | | 0% | |========================================| 100%# define colour palette
st_colors = c( wes_palette("Zissou1", n = 12, type = "continuous")[c(12,12,12,1,1,1)], "white" )
# plot tree with ancestral states as pies
zx=plotAncStatesPie(t=tre,
pie_colors=st_colors,
node_labels_size=0,
node_pie_size = 0.75,
#node_pie_nudge=0.5,
tip_pie_size=0.3,
tip_pie_nudge_x = 0.75,
node_pie_nudge_x = 0.5,
tip_labels=FALSE,
xlim=c(0,106),
state_transparency = 0.75,
layout="circular")
zx
Transition rates
Here’s where we test the hypotheses. First, we want to know if FPs are more often descendants from SD lineages than SM lineages.
# read posterior
Hmod1 <- read.table("../output/HiSSE/PolyHiSSE_MAP_coen_r1.model.log", header = T)
Hmod2 <- read.table("../output/HiSSE/PolyHiSSE_MAP_coen_r2.model.log", header = T)
# burn in the first 40,000 iterations and combine
rdat <- rbind(Hmod1[401:4400,c(1,42:47)],Hmod2[401:4400,c(1,42:47)])
# what's the fraction of the posterior in which SM -> FP > SD -> FP
sum(rdat$transition_rates.1. - rdat$transition_rates.6. > 0) / length(rdat$Iteration)## [1] 0.04475# What's the mean and 95% HPD interval of the difference
H1 <- data.frame(mean = mean(rdat$transition_rates.6. - rdat$transition_rates.1.),
lwr = HPDinterval(mcmc((rdat$transition_rates.6. - rdat$transition_rates.1.)))[1],
upr = HPDinterval(mcmc((rdat$transition_rates.6. - rdat$transition_rates.1.)))[2])
H1 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| mean | lwr | upr |
|---|---|---|
| 0.0043783 | -0.000519 | 0.0098235 |
We’ll plot mean rates for all transitions in Fig 2.
rdat %>% pivot_longer(cols = c(2:7), names_to = "Transition", values_to = "Rate") %>%
group_by(Transition) %>%
summarise(mean = mean(Rate))We plot histograms of the transitions rates towards FP for Fig. 2c.
# reshape to long
rdat_long <- pivot_longer(rdat, cols = c(2,7), names_to = "Transition", values_to = "Rate")
# get colour palette
pal <- wes_palette("Zissou1", n = 12, type = "continuous")
# plot rates
p1 <- ggplot(rdat_long,aes(x = Rate, fill = Transition)) +
theme_bw(base_size = 8) +
labs(x = "Transition rate", y = "Posterior frequency") +
geom_histogram(
binwidth = 0.0004,
linewidth = 0.25,
alpha = 0.75,
colour = "grey20",
position = "identity"
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
#xlim(0, 1.0) +
scale_fill_manual(
breaks = c("transition_rates.1.", "transition_rates.6."),
values = pal[c(1,7)],
labels = c(
"SM",
"SD"
),
name = "Ancestor of FP"
) +
theme_minimal(base_size = 7) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p1
Second, we want to know whether SD or SM is the most common descendant of FP.
# read posterior
# what's the fraction of the posterior in which FP -> SM > FP -> SM
sum(rdat$transition_rates.3. - rdat$transition_rates.4. > 0) / length(rdat$Iteration)## [1] 0.473125# What's the mean and 95% HPD interval of the difference
H2 <- data.frame(mean = mean(rdat$transition_rates.4. - rdat$transition_rates.3.),
lwr = HPDinterval(mcmc((rdat$transition_rates.4. - rdat$transition_rates.3.)))[1],
upr = HPDinterval(mcmc((rdat$transition_rates.4. - rdat$transition_rates.3.)))[2])
H2 %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| mean | lwr | upr |
|---|---|---|
| 0.0006824 | -0.0067679 | 0.0088242 |
We plot histograms of the transitions rates away FP for Fig. 2d
# reshape to long
rdat_long <- pivot_longer(rdat, cols = c(4,5), names_to = "Transition", values_to = "Rate")
# get colour palette
pal <- wes_palette("Zissou1", n = 12, type = "continuous")
# plot rates
p2 <- ggplot(rdat_long, aes(x = Rate, fill = Transition)) +
theme_bw(base_size = 8) +
labs(x = "Transition rate", y = "Posterior frequency") +
geom_histogram(
binwidth = 0.0004,
linewidth = 0.25,
alpha = 0.75,
colour = "grey20",
position = "identity"
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
#xlim(0, 1.0) +
scale_fill_manual(
breaks = c("transition_rates.3.", "transition_rates.4."),
values = pal[c(1,7)],
labels = c(
"SM",
"SD"
),
name = "Descendant of FP"
) +
theme_minimal(base_size = 7) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p2
Create multipanel figure
lay_matrix <- rbind(c(1,1),
c(1,1),
c(2,3))
Fig2 <- grid.arrange(zz, p1, p2, layout_matrix = lay_matrix)
#ggsave("../../Writing/Figures/Fig2_V3.pdf", plot = Fig2, device = "pdf", units = "mm", width = 180, height = 160)We plot histograms of the transitions rates between SD and SM for the Supporting Material
# reshape to long
rdat_long <- pivot_longer(rdat, cols = c(3,6), names_to = "Transition", values_to = "Rate")
# get colour palette
pal <- wes_palette("Zissou1", n = 12, type = "continuous")
# plot rates
p3 <- ggplot(rdat_long,aes(x = Rate, fill = Transition)) +
theme_bw(base_size = 8) +
labs(x = "Transition rate", y = "Posterior frequency") +
geom_histogram(
binwidth = 0.0004,
linewidth = 0.5,
alpha = 0.75,
colour = "white",
position = "identity"
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
#xlim(0, 1.0) +
scale_fill_manual(
breaks = c("transition_rates.2.", "transition_rates.5."),
values = c("grey80", "grey20"),
labels = c(
"SM to SD",
"SD to SM"
),
name = "Transition"
) +
theme_minimal(base_size = 7) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p3
Events
Stochastic character mapping allows us to estimate the number of events for each type of transition.
# it doesn't much matter which run we use
events <- read.table("../output/HiSSE/PolyHiSSE_MAP_coen_r1_events.tsv", header = T)
events <- read.table("../output/HiSSE/PolyHiSSE_MAP_coen_r2_events.tsv", header = T)
SDtoFP <- events %>%
filter(end_state == 1 | end_state == 4) %>%
filter(start_state == 2 | start_state == 5) %>%
group_by(iteration) %>% summarise(n = sum(end_state))
SMtoFP <- events %>%
filter(end_state == 1 | end_state == 4) %>%
filter(start_state == 0 | start_state == 3) %>%
group_by(iteration) %>% summarise(n = sum(end_state))
FPtoSD <- events %>%
filter(end_state == 2 | end_state == 5) %>%
filter(start_state == 1 | start_state == 5) %>%
group_by(iteration) %>% summarise(n = sum(end_state))
FPtoSM <- events %>%
filter(end_state == 2 | end_state == 5) %>%
filter(start_state == 1 | start_state == 4) %>%
group_by(iteration) %>% summarise(n = sum(end_state))
FP_table <- data.frame(FP_transition = rep(c("gain", "loss"), each = 2),
from_to = rep(c("SD", "SM"),2),
median = c(median(SDtoFP$n),
median(SMtoFP$n),
median(FPtoSD$n),
median(FPtoSM$n)),
lwr = c(HPDinterval(mcmc(SDtoFP$n))[1],
HPDinterval(mcmc(SMtoFP$n))[1],
HPDinterval(mcmc(FPtoSD$n))[1],
HPDinterval(mcmc(FPtoSM$n))[1]),
upr = c(HPDinterval(mcmc(SDtoFP$n))[2],
HPDinterval(mcmc(SMtoFP$n))[2],
HPDinterval(mcmc(FPtoSD$n))[2],
HPDinterval(mcmc(FPtoSM$n))[2]))
FP_table %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| FP_transition | from_to | median | lwr | upr |
|---|---|---|---|---|
| gain | SD | 42 | 4 | 66 |
| gain | SM | 3 | 1 | 12 |
| loss | SD | 89 | 2 | 315 |
| loss | SM | 68 | 2 | 167 |
Diversification dynamics
We start by plotting diversification rates.
Speciation
# reshape to long
rdat <- rbind(Hmod1[401:4400,c(1,c(28:33,5:10))],Hmod2[401:4400,c(1,c(28:33,5:10))])
colnames(rdat) <-
c(
"Iteration",
"speciation_SM_0",
"speciation_FP_0",
"speciation_SD_0",
"speciation_SM_1",
"speciation_FP_1",
"speciation_SD_1",
"extinction_SM_0",
"extinction_FP_0",
"extinction_SD_0",
"extinction_SM_1",
"extinction_FP_1",
"extinction_SD_1"
)
rdat_long <- pivot_longer(rdat, cols = c(2:13), names_to = "Event", values_to = "Rate") %>% separate_wider_delim(Event, delim = "_", names=c("Type", "State", "Hidden"))
pal <- wes_palette("Zissou1", 12, type = "continuous")
d1 <- ggplot(data = rdat_long[rdat_long$Type == "speciation",], aes(x = Rate, fill = State, alpha = Hidden)) +
geom_histogram(bins = 80, colour = "white", position = "identity", size = 0.2) +
scale_fill_manual(values = pal[c(12,7,1)], name = "Colour state") +
scale_alpha_manual(values = c(0.5,1), name = "Hidden state") +
labs(y = "Posterior frequency", title = "(A)") +
theme_bw(base_size = 8) +
theme(panel.grid = element_blank()) +
xlim(0,0.8) #+
#facet_grid(cols = vars(Hidden), scales = "free")
d1
spe_wide <- rdat %>%
mutate(diff_0vs1_0 = speciation_SM_0 - speciation_FP_0) %>%
mutate(diff_0vs2_0 = speciation_SM_0 - speciation_SD_0) %>%
mutate(diff_2vs1_0 = speciation_SD_0 - speciation_FP_0) %>%
mutate(diff_0vs1_1 = speciation_SM_1 - speciation_FP_1) %>%
mutate(diff_0vs2_1 = speciation_SM_1 - speciation_SD_1) %>%
mutate(diff_2vs1_1 = speciation_SD_1 - speciation_FP_1) %>%
dplyr::select(Iteration, starts_with("diff"))
n <- nrow(spe_wide)
spe_long <- pivot_longer(spe_wide, cols = c(2:7), names_to = "contrast", values_to = "diff") %>%
separate_wider_delim(cols = contrast, delim = "_", names = c("prefix", "states", "hidden") ) %>%
dplyr::select(!prefix)
spe_test <- spe_long %>%
group_by(states, hidden) %>%
summarise(Mean = mean(diff),
Lwr = HPDinterval(mcmc(diff))[1],
Upr = HPDinterval(mcmc(diff))[2],
PMCMC = sum(diff > 0)/n)
spe_test %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| states | hidden | Mean | Lwr | Upr | PMCMC |
|---|---|---|---|---|---|
| 0vs1 | 0 | -0.0649483 | -0.1156369 | -0.0208205 | 0.000000 |
| 0vs1 | 1 | -0.2669787 | -0.4855430 | -0.0683127 | 0.000000 |
| 0vs2 | 0 | -0.0203889 | -0.0355874 | -0.0043077 | 0.003750 |
| 0vs2 | 1 | -0.0902799 | -0.1909930 | -0.0046078 | 0.003750 |
| 2vs1 | 0 | -0.0445594 | -0.0920571 | -0.0005582 | 0.009375 |
| 2vs1 | 1 | -0.1766988 | -0.3754330 | -0.0081780 | 0.009375 |
Extinction
d2 <- ggplot(data = rdat_long[rdat_long$Type == "extinction",], aes(x = Rate, fill = State, alpha = Hidden)) +
geom_histogram(bins = 80, colour = "white", position = "identity", size = 0.2) +
scale_fill_manual(values = pal[c(12,7,1)], name = "Colour state") +
scale_alpha_manual(values = c(0.5,1), name = "Hidden state") +
labs(y = "Posterior frequency", title = "(B)") +
theme_bw(base_size = 8) +
theme(panel.grid = element_blank()) +
xlim(0,0.2) +
ylim(0,500)
#facet_grid(cols = vars(Hidden), scales = "free")
d2
ext_wide <- rdat %>%
mutate(diff_0vs1_0 = extinction_SM_0 - extinction_FP_0) %>%
mutate(diff_0vs2_0 = extinction_SM_0 - extinction_SD_0) %>%
mutate(diff_2vs1_0 = extinction_SD_0 - extinction_FP_0) %>%
mutate(diff_0vs1_1 = extinction_SM_1 - extinction_FP_1) %>%
mutate(diff_0vs2_1 = extinction_SM_1 - extinction_SD_1) %>%
mutate(diff_2vs1_1 = extinction_SD_1 - extinction_FP_1) %>%
dplyr::select(Iteration, starts_with("diff"))
n <- nrow(ext_wide)
ext_long <- pivot_longer(ext_wide, cols = c(2:7), names_to = "contrast", values_to = "diff") %>%
separate_wider_delim(cols = contrast, delim = "_", names = c("prefix", "states", "hidden") ) %>%
dplyr::select(!prefix)
ext_test <- ext_long %>%
group_by(states, hidden) %>%
summarise(Mean = mean(diff),
Lwr = HPDinterval(mcmc(diff))[1],
Upr = HPDinterval(mcmc(diff))[2],
PMCMC = sum(diff > 0)/n)
ext_test %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| states | hidden | Mean | Lwr | Upr | PMCMC |
|---|---|---|---|---|---|
| 0vs1 | 0 | -0.0454779 | -0.1017026 | 0.0055703 | 0.064500 |
| 0vs1 | 1 | -0.0939605 | -0.2633435 | 0.0246104 | 0.064500 |
| 0vs2 | 0 | -0.0001481 | -0.0126795 | 0.0097767 | 0.482125 |
| 0vs2 | 1 | -0.0036169 | -0.0296310 | 0.0264651 | 0.482125 |
| 2vs1 | 0 | -0.0453298 | -0.1040458 | 0.0070959 | 0.075750 |
| 2vs1 | 1 | -0.0903436 | -0.2912837 | 0.0218272 | 0.075750 |
Net diversification
div_dat <- rdat %>%
mutate(diversification_SM_0 = get("speciation_SM_0") - get("extinction_SM_0")) %>%
mutate(diversification_FP_0 = get("speciation_FP_0") - get("extinction_FP_0")) %>%
mutate(diversification_SD_0 = get("speciation_SD_0") - get("extinction_SD_0")) %>%
mutate(diversification_SM_1 = get("speciation_SM_1") - get("extinction_SM_1")) %>%
mutate(diversification_FP_1 = get("speciation_FP_1") - get("extinction_FP_1")) %>%
mutate(diversification_SD_1 = get("speciation_SD_1") - get("extinction_SD_1")) %>%
dplyr::select(c(1,14:19))
div_dat_long <- pivot_longer(div_dat, cols = c(2:7), names_to = "Event", values_to = "Rate") %>% separate_wider_delim(Event, delim = "_", names=c("Type", "State", "Hidden"))
d3 <- ggplot(data = div_dat_long, aes(x = Rate, fill = State, alpha = Hidden)) +
geom_histogram(bins = 80, colour = "white", position = "identity", size = 0.2) +
scale_fill_manual(values = pal[c(12,7,1)], name = "Colour state") +
scale_alpha_manual(values = c(0.5,1), name = "Hidden state") +
labs(y = "Posterior frequency", title = "(C)") +
theme_bw(base_size = 8) +
theme(panel.grid = element_blank()) +
xlim(0,0.6) #+
#facet_grid(cols = vars(Hidden), scales = "free")
d3
div_wide <- div_dat %>%
mutate(diff_0vs1_0 = diversification_SM_0 - diversification_FP_0) %>%
mutate(diff_0vs2_0 = diversification_SM_0 - diversification_SD_0) %>%
mutate(diff_2vs1_0 = diversification_SD_0 - diversification_FP_0) %>%
mutate(diff_0vs1_1 = diversification_SM_1 - diversification_FP_1) %>%
mutate(diff_0vs2_1 = diversification_SM_1 - diversification_SD_1) %>%
mutate(diff_2vs1_1 = diversification_SD_1 - diversification_FP_1) %>%
dplyr::select(Iteration, starts_with("diff"))
n <- nrow(div_wide)
div_long <- pivot_longer(div_wide, cols = c(2:7), names_to = "contrast", values_to = "diff") %>%
separate_wider_delim(cols = contrast, delim = "_", names = c("prefix", "states", "hidden") ) %>%
dplyr::select(!prefix)
div_test <- div_long %>%
group_by(states, hidden) %>%
summarise(Mean = mean(diff),
Lwr = HPDinterval(mcmc(diff))[1],
Upr = HPDinterval(mcmc(diff))[2],
PMCMC = sum(diff > 0)/n)
div_test %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| states | hidden | Mean | Lwr | Upr | PMCMC |
|---|---|---|---|---|---|
| 0vs1 | 0 | -0.0194704 | -0.0483134 | 0.0067590 | 0.067125 |
| 0vs1 | 1 | -0.1730182 | -0.3988528 | 0.0139063 | 0.040125 |
| 0vs2 | 0 | -0.0202408 | -0.0348651 | -0.0054547 | 0.004250 |
| 0vs2 | 1 | -0.0866630 | -0.1937440 | -0.0057165 | 0.011375 |
| 2vs1 | 0 | 0.0007704 | -0.0271121 | 0.0279462 | 0.511625 |
| 2vs1 | 1 | -0.0863552 | -0.2860556 | 0.0738765 | 0.151000 |
Finally combine results for all diversification rates for a table in the Supporting Material.
rate_type <- rep(c("speciation", "extinction", "net diversification", "turnover"), each = 6)
div_all <- bind_rows(list(spe_test,ext_test,div_test)) %>% bind_cols(rate_type) %>%
rename(Rate_type = ...7) %>% rename(Comparison = states) %>% rename(Hidden = hidden)
div_all %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)
#write.table(div_all, "../output/HiSSE/div_contrasts.tsv", quote = F, row.names = F, sep = "\t")And combine figures
Ecological and demographic predictors of FPs
We used six models to investigate the relationships between ecological factors, demographic factors and the phylogenetic distribution of FP. These analyses can be structured into 4 questions. For each question, we used a multinomial mixed effects model fitted using the package MCMCglmm (Hadfield 2010).
Are FPs more common in temperate regions and open areas?
See introduction in Willink et al. (2024) for the rationale of this hypothesis. Model 1 uses the literature data to address this question. The model chunk was run on UPPMAX.
Read and prepare the data
# read in tree
tre <- read.nexus(file = "../data/processed/HiSSE/Willink_MAP_coen.tree")
# read biome data from Willink et al
biodat <- read.csv("../data/processed/glmm/biome_dat.csv", sep = ",", header = T)[,-c(2:5)]
# make latitudinal range variable
biodat$Lat <- factor(paste0(biodat$Biome_1, biodat$Biome_2, biodat$Biome_3))
# revalue levels both warm and cold go to temperate and tropical stays as tropical
biodat$Lat <- dplyr::recode(biodat$Lat, "T" = "Trp", C = "Tmp", "W" = "Tmp", "WC" = "Tmp" )
biodat$Lat[biodat$Lat == "TW" | biodat$Lat == "TWC" | biodat$Lat == ""] <- NA
biodat <- biodat[,c("Taxon", "Lat")]
# number of ambiguous taxa
sum(is.na(biodat$Lat))## [1] 83biodat$Lat <- droplevels(biodat$Lat)
#levels(as.factor(biodat$Lat))
# read in opennes data
opendat<-read.csv("../data/processed/glmm/Eco_lit_dat.csv", sep = ",", header = T, na.strings = c("", "NA"), colClasses = rep("factor", 4))[,-c(2,4)]
# number of ambiguous taxa
sum(is.na(opendat$Hab))## [1] 106General data description
## [1] 356## [1] 62## [1] 384## how many taxa with ambiguous habitat openness
length(ecodat$Hab[ecodat$Hab == "C/O"]) - sum(is.na(ecodat$Hab))## [1] 50## data prep
# make ambiguous habitats NA
ecodat$Hab[ecodat$Hab == "C/O"] <- NA
ecodat$Hab <- droplevels(as.factor(ecodat$Hab))
#keep a version with all data
ecodat_full <- ecodat
# drop species missing either ecological variable
ecodat <- ecodat[complete.cases(ecodat[,2:3]) == T,]
#write.table(ecodat, "../data/processed/glmm/Eco_lit_dat.tsv", quote = F, row.names = F, sep = "\t")
# number of taxa left
nrow(ecodat)## [1] 297Correlated evolution
First I have re-format character data to nex and prune the full Coenagrionidae tree for each analysis.
# drop species missing Latitude
ecodat_lat <- ecodat_full[complete.cases(ecodat_full[,2]) == T,]
# re-code latitudinal ranges as 0 = Tropical and 1 = Temperate
Lat <- as.factor(ecodat_lat$Lat)
levels(Lat) <- c("1","0")
names(Lat) <- ecodat_lat$Taxon
# write habitat openness data to nex
write.nexus.data(Lat, file = "../data/processed/CorrEvol/Lat_full.nex", interleaved = T, format = "standard")
# re-code colour states as 0 = monomrphic and 1 = polymorphic
FP <- as.factor(ecodat_lat$State)
levels(FP) <- c("0","1", "0")
names(FP) <- ecodat_lat$Taxon
# write female-colour data for correlation with latitude analysis
write.nexus.data(FP, file = "../data/processed/CorrEvol/FP_lat.nex", interleaved = T, format = "standard")
# check taxon labels match
matching <- name.check(tre, ecodat_lat, ecodat_lat$Taxon)
# prune tre to match data
lattre <- drop.tip(tre, matching$tree_not_data, rooted = TRUE)
# fix rounding errors to make ultrametric
lattre <- nnls.tree(cophenetic(lattre), lattre, rooted = TRUE, method = "ultrametric")
# write tree for latitude data
write.tree(lattre, "../data/processed/CorrEvol/Lat_full.tree")
# drop species missing openness
ecodat_hab <- ecodat_full[complete.cases(ecodat_full[,3]) == T,]
# re-code habitat data to 0 = closed and 1 = open
Hab <- as.factor(ecodat_hab$Hab)
levels(Hab) <- c("0","1")
names(Hab) <- ecodat_hab$Taxon
# write habitat openness data to nex
write.nexus.data(Hab, file = "../data/processed/CorrEvol/Hab_full.nex", interleaved = T, format = "standard")
# re-code colour states as 0 = monomrphic and 1 = polymorphic
FP <- as.factor(ecodat_hab$State)
levels(FP) <- c("0","1", "0")
names(FP) <- ecodat_hab$Taxon
# write female-colour data for correlation with habitat openness analysis
write.nexus.data(FP, file = "../data/processed/CorrEvol/FP_hab.nex", interleaved = T, format = "standard")
# check taxon labels match
matching <- name.check(tre, ecodat_hab, ecodat_hab$Taxon)
# prune tre to match data
habtre <- drop.tip(tre, matching$tree_not_data, rooted = TRUE)
# fix rounding errors to make ultrametric
habtre <- nnls.tree(cophenetic(habtre), habtre, rooted = TRUE, method = "ultrametric")
# write tree for latitude data
write.tree(habtre, "../data/processed/CorrEvol/Hab_full.tree")Evaluate model convergence
check_lat <- checkConvergence(list_files = c("~/Dropbox/Doctorado/COEN_Comparative/FP_Comparative/output/CorrEvol/CorrEvol_LatLat_FP_corr_RJ_run_1.log",
"~/Dropbox/Doctorado/COEN_Comparative/FP_Comparative/output/CorrEvol/CorrEvol_LatLat_FP_corr_RJ_run_2.log" ))## Reading in log file 1
## Reading in log file 1
## [1] "Calculating burn-in"
## [1] "Analyzing continuous parameters"## ACHIEVED CONVERGENCE
##
## BURN-IN SET AT 0
##
##
## CONTINUOUS PARAMETERS WITH NO VARIANTION AND EXCLUDED FROM CONVERGENCE ASSESSMENT
## RUN 1
## prob_gain_B_indep
##
## RUN 2
## prob_gain_B_indep
##
## LOWEST CONTINUOUS PARAMETER ESS
## RUN 1 -> rate_loss_B_when_A1 2506.67
## RUN 2 -> rate_loss_B_when_A1 2988.27
## ## means ESS
## prob_gain_A_indep 0.691202199 7633.020
## prob_loss_A_indep 0.659835041 7919.023
## prob_loss_B_indep 0.630842289 6539.327
## rate_gain_A_when_B0 0.004015526 7516.652
## rate_gain_A_when_B1 0.003125682 7516.887
## rate_gain_B_when_A0 0.001540409 7112.365
## rate_gain_B_when_A1 0.011193774 7682.339
## rate_loss_A_when_B0 0.002899544 7902.383
## rate_loss_A_when_B1 0.002321353 8002.000
## rate_loss_B_when_A0 0.003414729 6612.013
## rate_loss_B_when_A1 0.004668638 5494.934
## rf.1. 0.395266870 7760.988
## rf.2. 0.199671486 7625.444
## rf.3. 0.201768071 7686.252
## rf.4. 0.203293578 7539.165
check_hab <- checkConvergence(list_files = c("~/Dropbox/Doctorado/COEN_Comparative/FP_Comparative/output/CorrEvol/CorrEvol_HabHab_FP_corr_RJ_run_1.log",
"~/Dropbox/Doctorado/COEN_Comparative/FP_Comparative/output/CorrEvol/CorrEvol_HabHab_FP_corr_RJ_run_2.log" ))## Reading in log file 1
## Reading in log file 1
## [1] "Calculating burn-in"
## [1] "Analyzing continuous parameters"## ACHIEVED CONVERGENCE
##
## BURN-IN SET AT 0
##
##
## LOWEST CONTINUOUS PARAMETER ESS
## RUN 1 -> prob_loss_B_indep 2424.9
## RUN 2 -> prob_loss_B_indep 2143.34
## ## means ESS
## prob_gain_A_indep 0.5363659085 7830.279
## prob_gain_B_indep 0.0003749063 8002.000
## prob_loss_A_indep 0.9883779055 8002.000
## prob_loss_B_indep 0.9015246188 4568.239
## rate_gain_A_when_B0 0.0018600867 7965.379
## rate_gain_A_when_B1 0.0014710061 7895.614
## rate_gain_B_when_A0 0.0002846961 7717.625
## rate_gain_B_when_A1 0.0054126687 7796.415
## rate_loss_A_when_B0 0.0066152132 8002.000
## rate_loss_A_when_B1 0.0065827915 8002.000
## rate_loss_B_when_A0 0.0067073294 6120.550
## rate_loss_B_when_A1 0.0071488302 7226.067
## rf.1. 0.1988879226 7926.879
## rf.2. 0.2019202443 7950.462
## rf.3. 0.3962121198 7928.134
## rf.4. 0.2029797043 7836.778
Test for independence
For each environmental prediction, I used the RJ mixture to estimate the probability that habitat evolution is dependent of female colour state and the probability that the evolution of FLCPs depends on the habitat occupied by a lineage.
# read posterior parameter estimates for correlation with latitude
lat_par <- read.table("../output/CorrEvol/CorrEvol_LatLat_FP_corr_RJ.log", header = T, sep = "\t")
n <- nrow(lat_par)
# compute probabilities of independent evolution
Lat_gain_indep <- lat_par %>% summarise(p = sum(prob_gain_A_indep)/n)
Lat_loss_indep <- lat_par %>% summarise(p = sum(prob_loss_A_indep)/n)
Fem_lat_gain_indep <- lat_par %>% summarise(p = sum(prob_gain_B_indep)/n)
Fem_lat_loss_indep <- lat_par %>% summarise(p = sum(prob_loss_B_indep)/n)
# read posterior parameter estimates for correlation with habitat
hab_par <- read.table("../output/CorrEvol/CorrEvol_HabHab_FP_corr_RJ.log", header = T, sep = "\t")
n <- nrow(hab_par)
# compute probabilities of independent evolution
Hab_gain_indep <- hab_par %>% summarise(p = sum(prob_gain_A_indep)/n)
Hab_loss_indep <- hab_par %>% summarise(p = sum(prob_loss_A_indep)/n)
Fem_hab_gain_indep <- hab_par %>% summarise(p = sum(prob_gain_B_indep)/n)
Fem_hab_loss_indep <- hab_par %>% summarise(p = sum(prob_loss_B_indep)/n)
# summarise probabilities in data frame
pp_table <- data.frame(model = rep(c("Temperate range", "Open habitat"), each = 4),
Transition = rep(c("Gain", "Gain", "Loss", "Loss"),2),
Trait = rep(c("Environmental \nfactor", "Female-colour \npolymorphism"), 4),
value = rbind(Lat_gain_indep, Fem_lat_gain_indep, Lat_loss_indep, Fem_lat_loss_indep,
Hab_gain_indep, Fem_hab_gain_indep, Hab_loss_indep, Fem_hab_loss_indep))
# calculate probability for a raw Bayes factor of 100 (decisive evidence)
BF1 <- c(100)
P1 = BF1/(1+BF1)
# calculate probability for a raw Bayes factor of 10 (strong evidence)
BF2 <- c(10)
P2 = BF2/(1+BF2)
# calculate probability for a raw Bayes factor of 3.2 (substantial evidence)
BF3 <- c(3.2)
P3 = BF3/(1+BF3)
# define colour palette
pal <- wes_palette("Zissou1", 12, type = "continuous")
# plot probabilities by character and transition type
PP <- ggplot(pp_table, aes(x = Transition, y = 1-p,)) +
geom_bar(aes(fill = Trait), color = "gray20", stat = "identity", width = 0.4,
position = position_dodge(width = 0.5), alpha = 0.6) +
ylim(0,1) +
labs (y = "Probability of correlated evolution")+
facet_grid(cols = vars(model)) +
theme_light(base_size = 8) +
theme(panel.grid = element_blank()) +
theme(legend.position = "top") +
#theme(legend.key.spacing.y = unit(0.2, 'cm')) +
scale_fill_manual(
values = pal[c(5, 9)], name = "Trait")+
geom_hline(yintercept=0.5, linetype="solid", linewidth = 0.3, color = "gray50") +
geom_hline(yintercept=P3, linetype="longdash", linewidth = 0.3, color = "gray50") +
geom_hline(yintercept=1-P3, linetype="longdash", linewidth = 0.3, color = "gray50") +
geom_hline(yintercept=P2, linetype="dashed", linewidth = 0.3, color = "gray50") +
geom_hline(yintercept=1-P2, linetype="dashed", linewidth = 0.3, color = "gray50") +
geom_hline(yintercept=P1, linetype="dotted", linewidth = 0.3, color = "gray50") +
geom_hline(yintercept=1-P1, linetype="dotted", linewidth = 0.3, color = "gray50")
PP
#ggsave(PP, filename = "../../Writing/Figures/Corr_test.pdf", device = "pdf", width = 90, height = 60, units = "mm", dpi = 600)Compute Bayes factors for transitions with some support of independent or correlated evolution.
# Latitudinal range model
BF_lat_gain = (1 - Fem_lat_gain_indep)/ (1 - (1- Fem_lat_gain_indep)) # correlated
BF_lat_gain## p
## 1 Inf# Habitat openness model
BF_FP_hab_gain = (1 - Fem_hab_gain_indep)/(1 - (1 - Fem_hab_gain_indep)) # correlated
BF_FP_hab_gain## p
## 1 2666.333## p
## 1 85.04301## p
## 1 9.154822Ancestral state reconstructions
Plot ancestral state reconstructions for the correlated evolution model between latitudinal range and FP.
CHARACTER_A <- "Lat"
CHARACTER_B <- "FP"
STATE_LABELS <- c("0" = "Tropical - FM", "1" = "Tropical - FP", "2" = "Temperate - FM", "3" = "Temperate - FP")
tree_file <- paste0("../output/CorrEvol/CorrEvol_LatLat_FP_ase_corr_RJ.tree")
# process the ancestral states
ase <- processAncStates(tree_file,
state_labels = STATE_LABELS)## | | | 0% | |========================================| 100%pal <- wes_palette("Zissou1", 12, type = "continuous")
p_lat <- plotAncStatesPie(t = ase,
tip_labels = TRUE,
tip_labels_size = 0.8,
node_pie_size = 0.8,
tip_pie_size = 0.3,
tip_labels_offset = 1.5,
tip_labels_italics = TRUE,
tip_pie_nudge_x = 0.5,
pie_colors = pal[c(1,1,1,1,10,12,7,5)],
colour ="gray20",
size = 0.2) +
theme(legend.position = c(0.12,0.92),
legend.title=element_text(size=5),
legend.text = element_text(size=5),
legend.key.size = unit(3,"mm") )
p_lat
#ggsave(plot = p_lat, filename = paste0("../output/CorrEvol/",CHARACTER_A,"_",CHARACTER_B,"_ASE_corr_RJ_Pie.pdf"), device = "pdf", dpi = 600, units = "mm", width = 90, height = 210)Plot ancestral state reconstructions for the correlated evolution model between habitat openness and FP.
CHARACTER_A <- "Hab"
CHARACTER_B <- "FP"
STATE_LABELS <- c("0" = "Closed - FM", "1" = "Closed - FP", "2" = "Open - FM", "3" = "Open - FP")
tree_file <- paste0("../output/CorrEvol/CorrEvol_HabHab_FP_ase_corr_RJ.tree")
# process the ancestral states
ase <- processAncStates(tree_file,
state_labels = STATE_LABELS)## | | | 0% | |========================================| 100%p_hab <- plotAncStatesPie(t = ase,
tip_labels = TRUE,
tip_labels_size = 0.8,
node_pie_size = 0.8,
tip_pie_size = 0.3,
tip_labels_offset = 1.5,
tip_labels_italics = TRUE,
tip_pie_nudge_x = 0.5,
pie_colors = pal[c(1,1,1,10,12,1,7,5)],
colour ="gray20",
size = 0.2) +
theme(legend.position = c(0.12,0.92),
legend.title=element_text(size=5),
legend.text = element_text(size=5),
legend.key.size = unit(3,"mm") )
p_hab
Quantify event sequences
The MRCA of Coenagrionidae likely was female-monomorphic and evolved in tropical-open habitats. For transitions to temperate and polymorphic states, I want to quantify the relative frequencies of environment first vs colour-state first transitions using stochastic character mapping. To do this, I’m adapting code from Landis et al. (2021).
source("./CorrEvol/triplet_util.R")
# file system
fp = "../"
base_fn = c("CorrEvol_LatLat_FP_corr_RJ")
out_fp = paste0(fp, "output/CorrEvol/")
history_fn = paste0(out_fp, base_fn,".events.tsv")
model_fn = paste0(out_fp, base_fn,".log")
# process data
d_raw = list()
d_param = list()
d_trip = list()
# number od states is 2x2 in correlated evoltion model
n_states = 4
# get datasets
f_burn = 0
thinby = 1
# read triplet data
d_raw = make_dat_triplet_raw(history_fn, n_states, f_burn, thinby )
# read parameters
d_param = read.table(model_fn, sep="\t", header=T)
# transition matrix - transitions that cannot occur are coded as 0
graph <- list(matrix(data = c(1,1,1,0,1,1,0,1,1,0,1,1,0,1,1,1),
nrow = 4, byrow = TRUE,
dimnames = list(c("1","2","3","4"), c("1","2","3","4"))))
# compute triplets
d_trip = make_triplet_df(d_raw, d_param, graph[[1]], subroot=T)
# put environment first (S2 == 3) vs female colour first (S2 == 2) triplets into a table
Lat_trip <- d_trip %>% filter(S1 == 1 & S3 == 4) %>%
group_by(iteration) %>%
summarise(Lat_first = sum(S2 == 3),
FP_first = sum(S2 == 2)) %>%
pivot_longer(cols = c(Lat_first, FP_first), names_to = "Sequence", values_to = "Count") %>%
group_by(Sequence) %>%
summarise(median = median(Count),
Lwr = HPDinterval(mcmc(Count))[1],
Upr = HPDinterval(mcmc(Count))[2])Lat_trip %>% kbl(booktabs =T, align= c(rep('l',1),rep('r', 3))) %>%
kable_styling(latex_options ="striped")| Sequence | median | Lwr | Upr |
|---|---|---|---|
| FP_first | 1 | 0 | 5 |
| Lat_first | 33 | 27 | 38 |
Box plot of posterior event series.
seq_d <- d_trip %>% filter(S1 == 1 & S3 == 4) %>%
group_by(iteration) %>%
summarise(Lat_first = sum(S2 == 3),
FP_first = sum(S2 == 2)) %>%
pivot_longer(cols = c(Lat_first, FP_first), names_to = "Sequence", values_to = "Count")
P_ev <- ggplot(seq_d, aes(x = Sequence, y=Count)) +
geom_boxplot(colour = "grey30") +
labs(y = "Transitions to FP and temperate ranges", x = "Event sequence") +
scale_x_discrete(labels = c("FP first", "Temperate first") ) +
theme_light(base_size = 8) +
theme(panel.grid = element_blank())
P_ev
##ggsave(P_ev, filename = "../../Writing/Figures/Triplet_counts.pdf", device = "pdf", width = 55, height = 60, units = "mm", dpi = 600)####Transition rates
Finally, plot the state-dependent transition rates from the correlated evolution models.
# specify the input files
Lat_rates <- read.table("../output/CorrEvol/CorrEvol_LatLat_FP_corr_RJ.log",sep = "\t", header =T)[2003:8002,]
Hab_rates <- read.table("../output/CorrEvol/CorrEvol_HabHab_FP_corr_RJ.log", sep = "\t", header =T)[2003:8002,]
# reshape data from latitudinal corr model
Lat_rates_long <- Lat_rates %>%
dplyr::select(c(1,10:17)) %>%
pivot_longer( , cols = c(2:9), names_to = "Transition", values_to = "Rate") %>%
separate_wider_delim(Transition, delim = "_when_" , names=c("Transition", "Condition")) %>%
separate_wider_delim(Transition, delim = "_" , names=c("empty", "Transition", "Trait")) %>%
dplyr::select(-empty)
# plot rates
Lat_gain <- ggplot(data = Lat_rates_long[Lat_rates_long$Transition == "gain",], aes(x = Rate, fill = interaction(Trait,Condition))) +
geom_histogram(bins = 50, colour = "grey20", position = "identity", linewidth = 0.2, alpha = 0.7,) +
scale_fill_manual(values = pal[c(12,10,7,5)], name = "Transition",
labels = c("FM tropical to \nFP tropical", "FM temperate to \nFP temperate ", "FM tropical to \nFM temperate", "FP tropical to \nFP temperate")) +
guides(fill = guide_legend(nrow = 4, title.position="top") ) +
labs(y = "Posterior frequency", title = "(C)") +
theme_bw(base_size = 8) +
theme(panel.grid = element_blank(), legend.position = "right", legend.key.spacing.y = unit(0.1, 'cm'))
#Lat_gain
#ggsave(Lat_gain, filename = "../../Writing/Figures/Corr_Lat_FP_gain.pdf", device = "pdf", width = 75, height = 80, units = "mm", dpi = 600)
# reshape data from habitat corr model
Hab_rates_long <- Hab_rates %>%
dplyr::select(c(1,10:17)) %>%
pivot_longer( , cols = c(2:9), names_to = "Transition", values_to = "Rate") %>%
separate_wider_delim(Transition, delim = "_when_" , names=c("Transition", "Condition")) %>%
separate_wider_delim(Transition, delim = "_" , names=c("empty", "Transition", "Trait")) %>%
dplyr::select(-empty)
# plot rates
Hab_gain <- ggplot(data = Hab_rates_long[Hab_rates_long$Transition == "gain",], aes(x = Rate, fill = interaction(Trait,Condition))) +
geom_histogram(bins = 50, colour = "grey20", position = "identity", linewidth = 0.2, alpha = 0.7) +
scale_fill_manual(values = pal[c(12,10,7,5)], name = "Transition",
labels = c("FM closed to \nFP closed", "FM open to \nFP open ", "FM closed to \nFM open", "FP closed to \nFP open")) +
labs(y = "Posterior frequency", title = "(D)") +
guides(fill = guide_legend(nrow = 4, title.position="top") ) +
theme_bw(base_size = 8) +
theme(panel.grid = element_blank(), legend.position = "right", legend.key.spacing.y = unit(0.1, 'cm'))
#Hab_gain
#ggsave(Hab_gain, filename = "../../Writing/Figures/Corr_Hab_FP_gain.pdf", device = "pdf", width =55, height = 80, units = "mm", dpi = 600)
Corr_evol <- grid.arrange(Lat_gain, Hab_gain, ncol = 1)
#ggsave(Corr_evol, filename = "../../Writing/Figures/Corr_FP_gain.pdf", device = "pdf", width = 85, height = 115, units = "mm", dpi = 600)Quantify differences in habitat model
n <- nrow(Hab_rates)
Corr_hab_summ <- Hab_rates %>%
mutate(diff_FMtoFP_HOvsHC = rate_gain_B_when_A1 - rate_gain_B_when_A0) %>%
mutate(diff_FPtoFM_HOvsHC = rate_loss_B_when_A1 - rate_loss_B_when_A0) %>%
mutate(diff_HCtoHO_FMvsFP = rate_gain_A_when_B0 - rate_gain_A_when_B1) %>%
mutate(diff_HOtoHC_FMvsFP = rate_loss_A_when_B0 - rate_loss_A_when_B1) %>%
dplyr::select(Iteration, starts_with("diff")) %>%
pivot_longer( starts_with("diff"), names_to = "Transition", values_to = "diff") %>%
group_by(Transition) %>%
summarise(Median = median(diff),
Lwr = HPDinterval(mcmc(diff))[1],
Upr = HPDinterval(mcmc(diff))[2],
PMCMC= sum(diff <= 0)/n) %>%
separate(Transition, into = c("diff","Transition", "Comparison"), sep = "_") %>% dplyr::select(!diff)
Corr_hab_summ %>% kbl(booktabs =T, align= c(rep('l',2),rep('r', 4))) %>%
kable_styling(latex_options ="striped")| Transition | Comparison | Median | Lwr | Upr | PMCMC |
|---|---|---|---|---|---|
| FMtoFP | HOvsHC | 0.0050971 | 0.0029543 | 0.0074009 | 0.0003333 |
| FPtoFM | HOvsHC | 0.0000000 | 0.0000000 | 0.0046154 | 0.9093333 |
| HCtoHO | FMvsFP | 0.0000000 | -0.0011606 | 0.0030213 | 0.6396667 |
| HOtoHC | FMvsFP | 0.0000000 | 0.0000000 | 0.0000000 | 0.9890000 |
Quantify differences in latitude model
n <- nrow(Lat_rates)
Corr_lat_summ <- Lat_rates %>%
mutate(diff_FMtoFP_TmpvsTrp = rate_gain_B_when_A1 - rate_gain_B_when_A0) %>%
mutate(diff_FPtoFM_Tmpvstrp = rate_loss_B_when_A1 - rate_loss_B_when_A0) %>%
mutate(diff_TrptoTmp_FMvsFP = rate_gain_A_when_B0 - rate_gain_A_when_B1) %>%
mutate(diff_TmptoTrp_FMvsFP = rate_loss_A_when_B0 - rate_loss_A_when_B1) %>%
dplyr::select(Iteration, starts_with("diff")) %>%
pivot_longer( starts_with("diff"), names_to = "Transition", values_to = "diff") %>%
group_by(Transition) %>%
summarise(Median = median(diff),
Lwr = HPDinterval(mcmc(diff))[1],
Upr = HPDinterval(mcmc(diff))[2],
PMCMC= sum(diff <= 0)/n) %>%
separate(Transition, into = c("diff","Transition", "Comparison"), sep = "_") %>% dplyr::select(!diff)
Corr_lat_summ %>% kbl(booktabs =T, align= c(rep('l',2),rep('r', 4))) %>%
kable_styling(latex_options ="striped")| Transition | Comparison | Median | Lwr | Upr | PMCMC |
|---|---|---|---|---|---|
| FMtoFP | TmpvsTrp | 0.009579 | 0.0054242 | 0.0139346 | 0.0000000 |
| FPtoFM | Tmpvstrp | 0.000000 | 0.0000000 | 0.0059061 | 0.6403333 |
| TmptoTrp | FMvsFP | 0.000000 | -0.0007243 | 0.0038024 | 0.7133333 |
| TrptoTmp | FMvsFP | 0.000000 | 0.0000000 | 0.0042250 | 0.7025000 |
Model 1
I ran this model on UPPMAX
module load R_packages/4.1.1
date
# Run Rscript in a clean R instance, output a logfile
Rscript --vanilla --verbose ./glmm/glmm_mm1.R > slurm-${SLURM_JOBID}.Rout 2>&1
# append logfile to this scripts logfile
cat slurm-${SLURM_JOBID}.Rout >> slurm-${SLURM_JOBID}.out
# remove Rout log
rm slurm-${SLURM_JOBID}.Rout
dateRun the model
# we need a matrix of phylogenetic distances to account for shared evolutionary history
inv.phylo <- inverseA(ecotre, nodes = "TIPS", scale = TRUE)
# residual covariance matrix 1/J*(I+J), I and J are identity and unit matrices respectively
IJ <- (1/3) * (diag(2) + matrix(1, 2, 2))
# kronecker prior for fixed effect close to being flat for the 2 way:
#(0 vs.1 , 1 vs. 2 and 0 vs. 2) marginal probabilities within each fixed effect
pr1 <- list(
B = list(mu = rep(0, 8), V = diag(8) * (1.7 + pi ^ 2 / 3)),
R = list(V = IJ, fix = 1),
G = list(G1 = list(
V = IJ,
nu = 1000,
alpha.mu = rep(0, 2),
alpha.V = diag(2)
))
)
# run the model without intercept to estimate the probability of each state in each combination of latitude and habitat
mm1 <-
MCMCglmm(
State ~ at.level(Hab, 'O'):at.level(Lat, 'Trp'):trait +
at.level(Hab, 'O'):at.level(Lat, 'Tmp'):trait +
at.level(Hab, 'C'):at.level(Lat, 'Trp'):trait +
at.level(Hab, 'C'):at.level(Lat, 'Tmp'):trait - 1,
random = ~ us(trait):Taxon,
rcov = ~ us(trait):units,
ginverse = list(Taxon = inv.phylo$Ainv),
family = "categorical",
data = ecodat,
prior = pr1,
pl = TRUE,
slice = TRUE,
nitt = 330000,
burnin = 3000,
thin = 100,
verbose = f
)
#write.table(mm1$VCV, file = "../output/glmm/mm1_VCV.tsv", quote = F, row.names = F, sep = "\t")
#write.table(mm1$Sol, file = "../output/glmm/mm1_Sol.tsv", quote = F, row.names = F, sep = "\t")
#write.table(mm1$Liab, file = "../output/glmm/mm1_Liab.tsv", quote = F, row.names = F, sep = "\t")Get all intercepts
mm1_Sol <- read.table(file = "../output/glmm/mm1_Sol.tsv", header = T, sep = "\t")
# rescale the intercepts as if the residual covariance matrix was zero
Delta <- cbind(c(-1, 1, 0), c(-1, 0, 1))
c2 <- (16 * sqrt(3)/(15 * pi))^2
D <- ginv(Delta %*% t(Delta)) %*% Delta
# reminder of the residual covariance structure
IJ <- (1/3) * (diag(2) + matrix(1, 2, 2))
# Tropical open habitat
Int1 <- t(apply(mm1_Sol[, 1:2], 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Trp_O <- (mcmc(exp(Int1)/rowSums(exp(Int1))))
# Temperate open landscape
Int2 <- t(apply(cbind(mm1_Sol[, 3:4]), 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Tmp_O <- (mcmc(exp(Int2)/rowSums(exp(Int2))))
# Tropical closed habitat
Int3 <- t(apply(cbind(mm1_Sol[, 5:6]), 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Trp_C <- (mcmc(exp(Int3)/rowSums(exp(Int3))))
# Temperate closed habitat
Int4 <- t(apply(cbind(mm1_Sol[, 7:8]), 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Tmp_C <- (mcmc(exp(Int4)/rowSums(exp(Int4))))Summary table
Summarise estimates in a table
Eco_FP_lit <- data.frame(Latitude =rep( c("tropical", "temperate"), each = 2),
Habitat = rep( c("open", "closed"), 2),
SM_mean = c(mean(Trp_O[,1]), mean(Trp_C[,1]), mean(Tmp_O[,1]), mean(Tmp_C[,1])),
FP_mean = c(mean(Trp_O[,2]), mean(Trp_C[,2]), mean(Tmp_O[,2]), mean(Tmp_C[,2])),
SD_mean = c(mean(Trp_O[,3]), mean(Trp_C[,3]), mean(Tmp_O[,3]), mean(Tmp_C[,3])))
Eco_FP_lit %>% kbl(booktabs =T, align= c(rep('l',2),rep('r', 3))) %>%
kable_styling(latex_options ="striped")| Latitude | Habitat | SM_mean | FP_mean | SD_mean |
|---|---|---|---|---|
| tropical | open | 0.3719972 | 0.0955807 | 0.5324222 |
| tropical | closed | 0.5640208 | 0.0391222 | 0.3968570 |
| temperate | open | 0.1136117 | 0.4535288 | 0.4328595 |
| temperate | closed | 0.3012055 | 0.3611345 | 0.3376600 |
Plot
We plot the results separately for each state. First prepare the data for plotting.
pal <- wes_palette("Zissou1", 12, type = "continuous")
n = nrow(mm1_Sol)
# female polymorphism
Plotdat<-data.frame(Lat = rep(c("Tropical", "Temperate", "Tropical", "Temperate"), each = n ),
Hab = rep(c("O", "C"), each = 2*n),
SM_p = c(Trp_O[,1], Tmp_O[,1], Trp_C[,1], Tmp_C[,1]),
FP_p = c(Trp_O[,2], Tmp_O[,2], Trp_C[,2], Tmp_C[,2]),
SD_p = c(Trp_O[,3], Tmp_O[,3], Trp_C[,3], Tmp_C[,3]))
Plotdat$Sample<-rownames(Plotdat)
Plotdat$Eco<-paste(Plotdat$Lat, Plotdat$Hab, sep="_")
#head(Plotdat)Ecological effects on the occurrence of sexual monomorphism (SM)
p1.0 <-
ggplot(Plotdat, aes(x = SM_p, color = Eco, fill = Eco)) + theme_bw(base_size = 8) +
labs(title = "(A)",
x = "Probability of sexual monomorphism", y = "Posterior frequency") +
geom_histogram(
binwidth = 0.02,
linewidth = 0.25,
alpha = 0.5,
colour = "grey20",
position = "identity",
aes(y = after_stat(count))
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
#xlim(0, 1.0) +
scale_fill_manual(
breaks = c("Tropical_C", "Tropical_O", "Temperate_C", "Temperate_O"),
values = pal[c(12, 8, 5, 1)],
labels = c(
"tropical-closed",
"tropical-open",
"temperate-closed",
"temperate-open"
),
name = "Latitude-habitat"
) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p1.0
Ecological effects on the occurrence of female polymorphism (FP)
p1.1 <-
ggplot(Plotdat, aes(x = FP_p, color = Eco, fill = Eco)) + theme_bw(base_size = 8) +
labs(title = "",
x = "Probability of female polymorphism", y = "Posterior frequency") +
geom_histogram(
binwidth = 0.02,
linewidth = 0.25,
alpha = 0.5,
colour = "grey20",
position = "identity",
aes(y = after_stat(count))
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
#xlim(0, 1.0) +
scale_fill_manual(
breaks = c("Tropical_C", "Tropical_O", "Temperate_C", "Temperate_O"),
values = pal[c(12, 8, 5, 1)],
labels = c(
"tropical-closed",
"tropical-open",
"temperate-closed",
"temperate-open"
),
name = "Latitude-habitat"
) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p1.1
#ggsave(p1.1,
# filename = "~/Dropbox/Doctorado/COEN_Comparative/Writing/Figures/Model1_FP.pdf", device = "pdf", units = "mm", width = 80, height = 40)Ecological effects on the occurrence of sexual dimorphism (SD)
p1.2 <-
ggplot(Plotdat, aes(x = SD_p, color = Eco, fill = Eco)) + theme_bw(base_size = 8) +
labs(title = "(B)",
x = "Probability of sexual dimorphism", y = "Posterior frequency") +
geom_histogram(
binwidth = 0.02,
linewidth = 0.25,
alpha = 0.5,
colour = "grey20",
position = "identity",
aes(y = after_stat(count))
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
#xlim(0, 1.0) +
scale_fill_manual(
breaks = c("Tropical_C", "Tropical_O", "Temperate_C", "Temperate_O"),
values = pal[c(12, 8, 5, 1)],
labels = c(
"tropical-closed",
"tropical-open",
"temperate-closed",
"temperate-open"
),
name = "Latitude-habitat"
) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p1.2
Export SM and SD for Supporting Material
Hypothesis test
FPs are more frequent in temperate regions, and in open areas of tropical regions. Sexual monomorphism is more common in closed, tropical habitats.
# length of posterior sample
n = length(Tmp_C[,2])
# Temperate vs tropical in open
LatO_m <- mean(Tmp_O[,2] - Trp_O[,2])
LatO_l <- HPDinterval(Tmp_O[,2] - Trp_O[,2])[1]
LatO_u <- HPDinterval(Tmp_O[,2] - Trp_O[,2])[2]
LatO_p <- sum(Trp_O[,2] - Tmp_O[,2] > 0) / n
# Temperate vs tropical in closed
LatC_m <- mean(Tmp_C[,2] - Trp_C[,2])
LatC_l <- HPDinterval(Tmp_C[,2] - Trp_C[,2])[1]
LatC_u <- HPDinterval(Tmp_C[,2] - Trp_C[,2])[2]
LatC_p <- sum(Trp_C[,2] - Tmp_C[,2] > 0) / n
# Open vs closed in tropics
HabTr_m <- mean(Trp_O[,2] - Trp_C[,2])
HabTr_l <- HPDinterval(Trp_O[,2] - Trp_C[,2])[1]
HabTr_u <- HPDinterval(Trp_O[,2] - Trp_C[,2])[2]
HabTr_p <- sum(Trp_C[,2] - Trp_O[,2] > 0) / n
# Open vs closed in tropics
HabTm_m <- mean(Tmp_O[,2] - Tmp_C[,2])
HabTm_l <- HPDinterval(Tmp_O[,2] - Tmp_C[,2])[1]
HabTm_u <- HPDinterval(Tmp_O[,2] - Tmp_C[,2])[2]
HabTm_p <- 1 - (sum(Tmp_O[,2] - Tmp_C[,2] > 0) / n)
# combine in dataframe
FP_tests <- data.frame(Comparsion = c("temperate-open vs tropical-open",
"temperate-closed vs tropical-closed",
"tropical-open vs tropical-closed",
"temperate-open vs temperate-closed"),
Mean = c(LatO_m, LatC_m, HabTr_m, HabTm_m),
Lower = c(LatO_l, LatC_l, HabTr_l, HabTm_l),
Upper = c(LatO_u, LatC_u, HabTr_u, HabTm_u),
PMCMC = c(LatO_p, LatC_p, HabTr_p, HabTm_p))
FP_tests %>%
kbl(booktabs =T, align= c(rep('l',1),rep('r', 4))) %>%
kable_styling(latex_options ="striped")| Comparsion | Mean | Lower | Upper | PMCMC |
|---|---|---|---|---|
| temperate-open vs tropical-open | 0.3579481 | 0.1151845 | 0.5922740 | 0.0003058 |
| temperate-closed vs tropical-closed | 0.3220123 | 0.0601452 | 0.6529775 | 0.0009174 |
| tropical-open vs tropical-closed | 0.0564585 | -0.0249451 | 0.1656533 | 0.0773700 |
| temperate-open vs temperate-closed | 0.0923943 | -0.2099826 | 0.4155171 | 0.2749235 |
# length of posterior sample
n = length(Tmp_C[,1])
# Temperate vs tropical in open
LatO_m <- mean(Tmp_O[,1] - Trp_O[,1])
LatO_l <- HPDinterval(Tmp_O[,1] - Trp_O[,1])[1]
LatO_u <- HPDinterval(Tmp_O[,1] - Trp_O[,1])[2]
LatO_p <- 1 - (sum(Trp_O[,1] - Tmp_O[,1] > 0) / n)
# Temperate vs tropical in closed
LatC_m <- mean(Tmp_C[,1] - Trp_C[,1])
LatC_l <- HPDinterval(Tmp_C[,1] - Trp_C[,1])[1]
LatC_u <- HPDinterval(Tmp_C[,1] - Trp_C[,1])[2]
LatC_p <- 1 - (sum(Trp_C[,1] - Tmp_C[,1] > 0) / n)
# Open vs closed in tropics
HabTr_m <- mean(Trp_O[,1] - Trp_C[,1])
HabTr_l <- HPDinterval(Trp_O[,1] - Trp_C[,1])[1]
HabTr_u <- HPDinterval(Trp_O[,1] - Trp_C[,1])[2]
HabTr_p <- 1 - (sum(Trp_C[,1] - Trp_O[,1] > 0) / n)
# Open vs closed in tropics
HabTm_m <- mean(Tmp_O[,1] - Tmp_C[,1])
HabTm_l <- HPDinterval(Tmp_O[,1] - Tmp_C[,1])[1]
HabTm_u <- HPDinterval(Tmp_O[,1] - Tmp_C[,1])[2]
HabTm_p <- sum(Tmp_O[,1] - Tmp_C[,1] > 0) / n
# combine in dataframe
SM_tests <- data.frame(Comparsion = c("temperate-open vs tropical-open",
"temperate-closed vs tropical-closed",
"tropical-open vs tropical-closed",
"temperate-open vs temperate-closed"),
Mean = c(LatO_m, LatC_m, HabTr_m, HabTm_m),
Lower = c(LatO_l, LatC_l, HabTr_l, HabTm_l),
Upper = c(LatO_u, LatC_u, HabTr_u, HabTm_u),
PMCMC = c(LatO_p, LatC_p, HabTr_p, HabTm_p))
SM_tests %>%
kbl(booktabs =T, align= c(rep('l',1),rep('r', 4))) %>%
kable_styling(latex_options ="striped")| Comparsion | Mean | Lower | Upper | PMCMC |
|---|---|---|---|---|
| temperate-open vs tropical-open | -0.2583855 | -0.4700423 | -0.0763370 | 0.0039755 |
| temperate-closed vs tropical-closed | -0.2628154 | -0.5421785 | 0.0462706 | 0.0489297 |
| tropical-open vs tropical-closed | -0.1920237 | -0.4099427 | -0.0076624 | 0.0327217 |
| temperate-open vs temperate-closed | -0.1875938 | -0.4551042 | 0.0344990 | 0.0529052 |
# length of posterior sample
n = length(Tmp_C[,3])
# Temperate vs tropical in open
LatO_m <- mean(Tmp_O[,3] - Trp_O[,3])
LatO_l <- HPDinterval(Tmp_O[,3] - Trp_O[,3])[1]
LatO_u <- HPDinterval(Tmp_O[,3] - Trp_O[,3])[2]
LatO_p <- 1 - (sum(Trp_O[,3] - Tmp_O[,3] > 0) / n)
# Temperate vs tropical in closed
LatC_m <- mean(Tmp_C[,3] - Trp_C[,3])
LatC_l <- HPDinterval(Tmp_C[,3] - Trp_C[,3])[1]
LatC_u <- HPDinterval(Tmp_C[,3] - Trp_C[,3])[2]
LatC_p <- 1 - (sum(Trp_C[,3] - Tmp_C[,3] > 0) / n)
# Open vs closed in tropics
HabTr_m <- mean(Trp_O[,3] - Trp_C[,3])
HabTr_l <- HPDinterval(Trp_O[,3] - Trp_C[,3])[1]
HabTr_u <- HPDinterval(Trp_O[,3] - Trp_C[,3])[2]
HabTr_p <- sum(Trp_C[,3] - Trp_O[,3] > 0) / n
# Open vs closed in tropics
HabTm_m <- mean(Tmp_O[,3] - Tmp_C[,3])
HabTm_l <- HPDinterval(Tmp_O[,3] - Tmp_C[,3])[1]
HabTm_u <- HPDinterval(Tmp_O[,3] - Tmp_C[,3])[2]
HabTm_p <- 1 - (sum(Tmp_O[,3] - Tmp_C[,3] > 0) / n)
# combine in dataframe
SD_tests <- data.frame(Comparsion = c("temperate-open vs tropical-open",
"temperate-closed vs tropical-closed",
"tropical-open vs tropical-closed",
"temperate-open vs temperate-closed"),
Mean = c(LatO_m, LatC_m, HabTr_m, HabTm_m),
Lower = c(LatO_l, LatC_l, HabTr_l, HabTm_l),
Upper = c(LatO_u, LatC_u, HabTr_u, HabTm_u),
PMCMC = c(LatO_p, LatC_p, HabTr_p, HabTm_p))
SD_tests %>%
kbl(booktabs =T, align= c(rep('l',1),rep('r', 4))) %>%
kable_styling(latex_options ="striped")| Comparsion | Mean | Lower | Upper | PMCMC |
|---|---|---|---|---|
| temperate-open vs tropical-open | -0.0995627 | -0.3614189 | 0.1427201 | 0.2293578 |
| temperate-closed vs tropical-closed | -0.0591969 | -0.3809587 | 0.2892279 | 0.3422018 |
| tropical-open vs tropical-closed | 0.1355652 | -0.0870159 | 0.3455783 | 0.1091743 |
| temperate-open vs temperate-closed | 0.0951994 | -0.2365086 | 0.4180145 | 0.2773700 |
Combine tables for Supporting Material
Repeat with field data only (Model 2)
In Model 2, we make sure this general pattern holds for the subset of data collected in the field.
Read and prepare the data
# Read in lit data
dat<-read.csv("../data/processed/glmm/Field_comp_dat.csv", sep = ",", header = T, na.strings = c("", "NA"), colClasses = c(rep("factor", 8), rep("numeric",3)))
# drop species missing either ecological variable
dat <- dat[complete.cases(dat[,6:8]) == T,]
# number of taxa left
nrow(dat)## [1] 244# prune data not in tree (featherlegs)
matching <- name.check(tre, dat, dat$Taxon)
ecodat2 <- dat[!(dat$Taxon %in% matching$data_not_tree),]
# prune trees to match data
ecotre2<-drop.tip(tre, matching$tree_not_data, rooted = TRUE)
# fix rounding errors to make ultrametric
ecotre2 <- nnls.tree(cophenetic(ecotre2), ecotre2, rooted = TRUE, method = "ultrametric")Run the model
# we need a matrix of phylogenetic distances to account for shared evolutionary history
inv.phylo <- inverseA(ecotre2, nodes = "TIPS", scale = TRUE)
# residual covariance matrix 1/J*(I+J), I and J are identity and unit matrices respectively
IJ <- (1/3) * (diag(2) + matrix(1, 2, 2))
# kronecker prior for fixed effect close to being flat for the 2 way:
#(0 vs.1 , 1 vs. 2 and 0 vs. 2) marginal probabilities within each fixed effect
pr1 <- list(
B = list(mu = rep(0, 8), V = diag(8) * (1.7 + pi ^ 2 / 3)),
R = list(V = IJ, fix = 1),
G = list(
G1 = list(
V = IJ,
nu = 1000,
alpha.mu = rep(0, 2),
alpha.V = diag(2)
),
G2 = list(
V = 1,
nu = 1000,
alpha.mu = rep(0, 1),
alpha.V = diag(1)
)
)
)
# run the model without intercept to estimate the probability of each state in each combination of latitude and habitat
# with SM as baseline for log-odds
mm2.1 <-
MCMCglmm(
relevel(FemState, ref="0") ~ at.level(Hab, 'O'):at.level(Lat, 'Trp'):trait +
at.level(Hab, 'O'):at.level(Lat, 'Tmp'):trait +
at.level(Hab, 'C'):at.level(Lat, 'Trp'):trait +
at.level(Hab, 'C'):at.level(Lat, 'Tmp'):trait - 1,
random = ~ us(trait):Taxon + Site,
rcov = ~ us(trait):units,
ginverse = list(Taxon = inv.phylo$Ainv),
family = "categorical",
data = ecodat2,
prior = pr1,
pl = TRUE,
slice = TRUE,
nitt = 330000,
burnin = 30000,
thin = 100,
verbose = F
)
# with SD as baseline for log-odds
mm2.2 <-
MCMCglmm(
relevel(FemState, ref="2") ~ at.level(Hab, 'O'):at.level(Lat, 'Trp'):trait +
at.level(Hab, 'O'):at.level(Lat, 'Tmp'):trait +
at.level(Hab, 'C'):at.level(Lat, 'Trp'):trait +
at.level(Hab, 'C'):at.level(Lat, 'Tmp'):trait - 1,
random = ~ us(trait):Taxon + Site,
rcov = ~ us(trait):units,
ginverse = list(Taxon = inv.phylo$Ainv),
family = "categorical",
data = ecodat2,
prior = pr1,
pl = TRUE,
slice = TRUE,
nitt = 330000,
burnin = 30000,
thin = 100,
verbose = F
)Get all intercepts
# rescale the intercepts as if the residual covariance matrix was zero
Delta <- cbind(c(-1, 1, 0), c(-1, 0, 1))
c2 <- (16 * sqrt(3)/(15 * pi))^2
D <- ginv(Delta %*% t(Delta)) %*% Delta
# Tropical open habitat
Int1 <- t(apply(mm2.1$Sol[, 1:2], 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Trp_O<- (mcmc(exp(Int1)/rowSums(exp(Int1))))
# Temperate open landscape
Int2 <- t(apply(cbind(mm2.1$Sol[, 3:4]), 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Tmp_O<-(mcmc(exp(Int2)/rowSums(exp(Int2)))) #
# Tropical closed habitat
Int3 <- t(apply(cbind(mm2.1$Sol[, 5:6]), 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Trp_C<-(mcmc(exp(Int3)/rowSums(exp(Int3))))
# Temperate closed habitat
Int4 <- t(apply(cbind(mm2.1$Sol[, 7:8]), 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Tmp_C<-(mcmc(exp(Int4)/rowSums(exp(Int4))))Summary table
Eco_FP_field <- data.frame(Latitude =rep( c("Tropical", "Temperate"), each = 2),
Habitat = rep( c("open", "closed"), 2),
SM_mean = c(mean(Trp_O[,1]), mean(Trp_C[,1]), mean(Tmp_O[,1]), mean(Tmp_C[,1])),
SD_mean = c(mean(Trp_O[,3]), mean(Trp_C[,3]), mean(Tmp_O[,3]), mean(Tmp_C[,3])),
FP_mean = c(mean(Trp_O[,2]), mean(Trp_C[,2]), mean(Tmp_O[,2]), mean(Tmp_C[,2])))
Eco_FP_field %>% kbl(booktabs =T, align= c(rep('l',2),rep('r', 3))) %>%
kable_styling(latex_options ="striped")| Latitude | Habitat | SM_mean | SD_mean | FP_mean |
|---|---|---|---|---|
| Tropical | open | 0.1053241 | 0.7247662 | 0.1699097 |
| Tropical | closed | 0.5821299 | 0.3104929 | 0.1073772 |
| Temperate | open | 0.0377187 | 0.1251959 | 0.8370855 |
| Temperate | closed | 0.5032556 | 0.3883073 | 0.1084371 |
Eco_FP <- rbind(Eco_FP_lit, Eco_FP_field)
Eco_FP$Source <- rep (c("literature", "field"), each = 4)
Eco_FP <- Eco_FP[,c(6,1:5)]
Eco_FP %>% kbl(booktabs =T, align= c(rep('l',3),rep('r', 3))) %>%
kable_styling(latex_options ="striped")| Source | Latitude | Habitat | SM_mean | FP_mean | SD_mean |
|---|---|---|---|---|---|
| literature | tropical | open | 0.3719972 | 0.0955807 | 0.5324222 |
| literature | tropical | closed | 0.5640208 | 0.0391222 | 0.3968570 |
| literature | temperate | open | 0.1136117 | 0.4535288 | 0.4328595 |
| literature | temperate | closed | 0.3012055 | 0.3611345 | 0.3376600 |
| field | Tropical | open | 0.1053241 | 0.1699097 | 0.7247662 |
| field | Tropical | closed | 0.5821299 | 0.1073772 | 0.3104929 |
| field | Temperate | open | 0.0377187 | 0.8370855 | 0.1251959 |
| field | Temperate | closed | 0.5032556 | 0.1084371 | 0.3883073 |
Plot
We plot the results separately for each state. First prepare the data for plotting.
pal <- wes_palette("Zissou1", 12, type = "continuous")
n = nrow(mm2.1$Sol)
#Female polymorphism
Plotdat<-data.frame(Lat = rep(c("Tropical", "Temperate", "Tropical", "Temperate"), each = n ),
Hab = rep(c("O", "C"), each = 2*n),
SM_p = c(Trp_O[,1], Tmp_O[,1], Trp_C[,1], Tmp_C[,1]),
FP_p = c(Trp_O[,2], Tmp_O[,2], Trp_C[,2], Tmp_C[,2]),
SD_p = c(Trp_O[,3], Tmp_O[,3], Trp_C[,3], Tmp_C[,3]))
Plotdat$Sample<-rownames(Plotdat)
Plotdat$Eco<-paste(Plotdat$Lat, Plotdat$Hab, sep="_")p2.0 <-
ggplot(Plotdat, aes(x = SM_p, color = Eco, fill = Eco)) + theme_bw(base_size = 8) +
labs(title = "",
x = "Probability of sexual monomorphism", y = "Posterior frequency") +
geom_histogram(
binwidth = 0.02,
linewidth = 0.25,
alpha = 0.5,
colour = "grey20",
position = "identity",
aes(y = after_stat(count))
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
#xlim(0, 1.0) +
scale_fill_manual(
breaks = c("Tropical_C", "Tropical_O", "Temperate_C", "Temperate_O"),
values = pal[c(12, 8, 5, 1)],
labels = c(
"Tropical-closed",
"Tropical-open",
"Temperate-closed",
"Temperate-open"
),
name = "Distribution-habitat"
) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p2.0
p2.1 <-
ggplot(Plotdat, aes(x = FP_p, color = Eco, fill = Eco)) + theme_bw(base_size = 8) +
labs(title = "",
x = "Probability of female polymorphism", y = "Posterior frequency") +
geom_histogram(
binwidth = 0.02,
linewidth = 0.25,
alpha = 0.5,
colour = "grey20",
position = "identity",
aes(y = after_stat(count))
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
#xlim(0, 1.0) +
scale_fill_manual(
breaks = c("Tropical_C", "Tropical_O", "Temperate_C", "Temperate_O"),
values = pal[c(12, 8, 5, 1)],
labels = c(
"Tropical-closed",
"Tropical-open",
"Temperate-closed",
"Temperate-open"
),
name = "Latitude-habitat"
) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p2.1
#ggsave(p2.1,
# filename = "~/Dropbox/Doctorado/COEN_Comparative/Writing/Figures/Model2_FP.pdf", device = "pdf", units = "mm", width = 80, height = 40)p2.2 <-
ggplot(Plotdat, aes(x = SD_p, color = Eco, fill = Eco)) + theme_bw(base_size = 8) +
labs(title = "",
x = "Probability of sexual dimorphism", y = "Posterior frequency") +
geom_histogram(
binwidth = 0.02,
linewidth = 0.25,
alpha = 0.5,
colour = "grey20",
position = "identity",
aes(y = after_stat(count))
) +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
#xlim(0, 1.0) +
scale_fill_manual(
breaks = c("Tropical_C", "Tropical_O", "Temperate_C", "Temperate_O"),
values = pal[c(12, 8, 5, 1)],
labels = c(
"Tropical-closed",
"Tropical-open",
"Temperate-closed",
"Temperate-open"
),
name = "Distribution-habitat"
) +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p2.2
Hypothesis test
# length of posterior sample
n = length(Tmp_C[,2])
# Temperate vs tropical in open
LatO_m <- mean(Tmp_O[,2] - Trp_O[,2])
LatO_l <- HPDinterval(Tmp_O[,2] - Trp_O[,2])[1]
LatO_u <- HPDinterval(Tmp_O[,2] - Trp_O[,2])[2]
LatO_p <- sum(Trp_O[,2] - Tmp_O[,2] > 0) / n
# Temperate vs tropical in closed
LatC_m <- mean(Tmp_C[,2] - Trp_C[,2])
LatC_l <- HPDinterval(Tmp_C[,2] - Trp_C[,2])[1]
LatC_u <- HPDinterval(Tmp_C[,2] - Trp_C[,2])[2]
LatC_p <- sum(Trp_C[,2] - Tmp_C[,2] > 0) / n
# Open vs closed in tropics
HabTr_m <- mean(Trp_O[,2] - Trp_C[,2])
HabTr_l <- HPDinterval(Trp_O[,2] - Trp_C[,2])[1]
HabTr_u <- HPDinterval(Trp_O[,2] - Trp_C[,2])[2]
HabTr_p <- sum(Trp_C[,2] - Trp_O[,2] > 0) / n
# Open vs closed in tropics
HabTm_m <- mean(Tmp_C[,1] - Tmp_O[,1])
HabTm_l <- HPDinterval(Tmp_C[,1] - Tmp_O[,1])[1]
HabTm_u <- HPDinterval(Tmp_C[,1] - Tmp_O[,1])[2]
HabTm_p <- sum(Tmp_O[,1] - Tmp_C[,1] > 0) / n
# combine in dataframe
FP_tests <- data.frame(comparsion = c("Temperate-open vs Tropical-open",
"Temperate-closed vs Tropical-closed",
"Tropical-open vs Tropical-closed",
"Temperate-open vs Temperate-closed"),
mean = c(LatO_m, LatC_m, HabTr_m, HabTm_m),
lower = c(LatO_l, LatC_l, HabTr_l, HabTm_l),
upper = c(LatO_u, LatC_u, HabTr_u, HabTm_u),
pMCMC = c(LatO_p, LatC_p, HabTr_p, HabTm_p))
FP_tests %>%
kbl(booktabs =T, align= c(rep('l',1),rep('r', 4))) %>%
kable_styling(latex_options ="striped")| comparsion | mean | lower | upper | pMCMC |
|---|---|---|---|---|
| Temperate-open vs Tropical-open | 0.6671758 | 0.3537082 | 0.9151402 | 0.000 |
| Temperate-closed vs Tropical-closed | 0.0010599 | -0.2758916 | 0.2829115 | 0.560 |
| Tropical-open vs Tropical-closed | 0.0625325 | -0.1570115 | 0.3484766 | 0.305 |
| Temperate-open vs Temperate-closed | 0.4655370 | 0.0270757 | 0.8434607 | 0.008 |
Are FLCPs more common in in species with high breeding density and male-biased OSR?
See introduction in Willink et al. (2024) for the rationale of this hypothesis.
Summary of field sampling
# How many breeding sites per country?
site_count <- aggregate(ecodat2$Site, by = list(ecodat2$Country), FUN = length)
site_count %>%
kbl() %>%
kable_material(c("striped", "hover"), full_width = F)| Group.1 | x |
|---|---|
| Argentina | 19 |
| Cameroon | 41 |
| Costa Rica | 48 |
| Cyprus | 17 |
| Guyana | 20 |
| NewZealand | 1 |
| Sweden | 37 |
| USA | 25 |
# Effort per site
site_time <- aggregate(ecodat2$Time, by = list(ecodat2$Site), FUN = sum)
effort <- data.frame(stat = c("min", "max", "mean", "sd"),
value = c(min(site_time$x/60),
max(site_time$x/60),
mean(site_time$x/60),
sd(site_time$x/60)))
effort %>%
kbl(booktabs =T, align= c('l','r')) %>%
kable_styling(latex_options ="striped")| stat | value |
|---|---|
| min | 0.200000 |
| max | 16.333333 |
| mean | 2.572917 |
| sd | 3.180338 |
Prepare the density data
How are the raw density and OSR data distributed
For Fig. 3c
pdat <- ggplot(data = ecodat2, aes(x = log(OSR), y = log(dens), colour = FemState)) +
geom_point(alpha = 0.7, size = 1) +
labs(x = "Log operational sex ratio", y = "Log adult density") +
guides(fill = guide_legend(ncol = 1, byrow = TRUE), color = NULL) +
scale_fill_manual(
breaks = c(0,1,2),
values = pal[c(1,12,7)],
labels = c("SM", "FP", "SD"),
name = "Colour state"
) +
scale_colour_manual(
breaks = c(0,1,2),
values = pal[c(1,12,7)],
labels = c("SM", "FP", "SD"),
name = "Colour state"
) +
theme_minimal(base_size = 8)
pdat 
Effects of adult density on the occurrence of FPs
In Model 3 we ask if the probability of FP increases with adult density at breeding sites.
Run the model
inv.phylo<-inverseA(ecotre2,nodes="TIPS",scale=TRUE)
# residual covariance matrix 1/J*(I+J), I and J are identity and unit matrices respectively
IJ <- (1/3) * (diag(2) + matrix(1, 2, 2))
# kronecker prior for fixed effect close to being flat for the 2 way:
#(0 vs.1 , 1 vs. 2 and 0 vs. 2) marginal probabilities within each fixed effect
prDem1 = list(
B = list(mu = rep(0, 4), V = diag(4) * (1.7 + pi ^ 2 / 3)),
R = list(V = IJ, fix = 1),
G = list(
G1 = list(
V = 1,
nu = 1000,
alpha.mu = rep(0, 1),
alpha.V = diag(1)
),
G2 = list(
V = 1,
nu = 1000,
alpha.mu = rep(0, 1),
alpha.V = diag(1)
)
)
)
# the model for the effect of density
# with SM as baseline
mm3.1 <- MCMCglmm(
relevel(FemState, ref = "0") ~ trait + log(dens):trait,
random = ~ Site + Taxon,
rcov = ~ us(trait):units,
ginverse = list(Taxon = inv.phylo$Ainv),
family = "categorical",
data = ecodat2,
prior = prDem1,
pl = TRUE,
slice = TRUE,
nitt = 330000,
burnin = 30000,
thin = 100,
verbose = F
)
# with SD as baseline
mm3.2 <- MCMCglmm(
relevel(FemState, ref = "2") ~ trait + log(dens):trait,
random = ~ Site + Taxon,
rcov = ~ us(trait):units,
ginverse = list(Taxon = inv.phylo$Ainv),
family = "categorical",
data = ecodat2,
prior = prDem1,
pl = TRUE,
slice = TRUE,
nitt = 330000,
burnin = 30000,
thin = 100,
verbose = F
)Get estimates in response scale
# rescale estimates assuming residual covariance matrix of zero
Delta <- cbind(c(-1, 1, 0), c(-1, 0, 1))
c2 <- (16 * sqrt(3)/(15 * pi))^2
D <- ginv(Delta %*% t(Delta)) %*% Delta
dens<-log(ecodat2$dens)
# make the estimate data frame
Pred<-data.frame(SM.mean = numeric(length = length(dens)),
FP.mean = numeric(length = length(dens)),
SD.mean = numeric(length = length(dens)),
SM.lwr = numeric(length = length(dens)),
FP.lwr = numeric(length = length(dens)),
SD.lwr = numeric(length = length(dens)),
SM.upr = numeric(length = length(dens)),
FP.upr = numeric(length = length(dens)),
SD.upr = numeric(length = length(dens)))
Pred$dens <- dens
Int<-data.frame(SM = numeric(length = 2*nrow(mm3.1$Sol)),
FP = numeric(length = 2*nrow(mm3.1$Sol)),
SD = numeric(length = 2*nrow(mm3.1$Sol)))
for (i in 1:length(dens)){
# rescale the estimates as if the residual covariance matrix was zero
Int1 <- t(apply((mm3.1$Sol[, 1:2]+dens[i]*mm3.1$Sol[, 3:4]), 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Int2 <- t(apply((mm3.2$Sol[, 1:2]+dens[i]*mm3.2$Sol[, 3:4]), 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Int <- rbind(Int1, Int2[,c(2,3,1)])
# median probability
Pred$SM.mean[i] <- mean(mcmc(exp(Int[,1])/rowSums(exp(Int))))
Pred$FP.mean[i] <- mean(mcmc(exp(Int[,2])/rowSums(exp(Int))))
Pred$SD.mean[i] <- mean(mcmc(exp(Int[,3])/rowSums(exp(Int))))
# lower HPD interval
Pred$SM.lwr[i] <- HPDinterval(mcmc(exp(Int[,1])/rowSums(exp(Int))))[1]
Pred$FP.lwr[i] <- HPDinterval(mcmc(exp(Int[,2])/rowSums(exp(Int))))[1]
Pred$SD.lwr[i] <- HPDinterval(mcmc(exp(Int[,3])/rowSums(exp(Int))))[1]
# upper HPD interval
Pred$SM.upr[i] <- HPDinterval(mcmc(exp(Int[,1])/rowSums(exp(Int))))[2]
Pred$FP.upr[i] <- HPDinterval(mcmc(exp(Int[,2])/rowSums(exp(Int))))[2]
Pred$SD.upr[i] <- HPDinterval(mcmc(exp(Int[,3])/rowSums(exp(Int))))[2]
}Plot predicted
# melt for plotting
Pred$sample <- rownames(Pred)
Pdat <- data.frame(dens = rep(dens,3))
Pdat$FemState <- rep(c("SM","FP", "SD"), each = length(dens))
Pdat$mean <- reshape::melt(Pred[,c(11,1:3)], id.vars = "sample")[,3]
Pdat$lwr <- reshape::melt(Pred[,c(11,4:6)], id.vars = "sample")[,3]
Pdat$upr <- reshape::melt(Pred[,c(11,7:9)], id.vars = "sample")[,3]
#head(Pdat)
# get species data points
Mdata <- data.frame(taxon = ecodat2$Taxon,
dens = log(ecodat2$dens),
SM = numeric(length(ecodat2$Taxon)),
FP = numeric(length(ecodat2$Taxon)),
SD = numeric(length(ecodat2$Taxon)))
for (i in 1:nrow(Mdata)){
if(ecodat2$FemState[i] == 0){
Mdata$SM[i] <- 1}
else{
if(ecodat2$FemState[i] == 1){
Mdata$FP[i] <-1}
else{
Mdata$SD[i] <-1
}
}
}
#head(Mdata)
# melt for plotting
MdataPlot <- reshape::melt(Mdata, id.vars = c("taxon", "dens"))
colnames(MdataPlot)[3:4] <- c("FemState", "mean")
pal <- wes_palette("Zissou1", 12, type = "continuous")
p3.2 <- ggplot(data = Pdat, aes(x=dens, y =mean, color = FemState ))+
geom_line(linewidth=0.5, show.legend = NA)+
theme_bw(base_size = 8)+
geom_ribbon(aes(ymax = upr, ymin = lwr, x=dens, fill=FemState),color=NA, alpha = 0.1)+
geom_point(data = MdataPlot, aes(color = FemState, fill = FemState), alpha = 0.3 , size = 1, position = position_jitter(height = 0,) )+
ylab(label = "Probability of \nsex-related colour state")+
xlab(label="Log adult density")+ # xlim(0,80)+
scale_color_manual(breaks = c("FP", "SD", "SM"),
values = c(pal[c(11,7,1)]),
name = "Colour state") +
scale_fill_manual(breaks = c("FP", "SD", "SM"),
values = c(pal[c(11,7,1)]),
name = "Colour state") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p3.2
Get log-odds estimates
dens<-log(ecodat2$dens)
# make the estimate data frame
Pred<-data.frame(FPvSM.mean = numeric(length = length(dens)),
FPvSD.mean = numeric(length = length(dens)),
FPvSM.lwr = numeric(length = length(dens)),
FPvSD.lwr = numeric(length = length(dens)),
FPvSM.upr = numeric(length = length(dens)),
FPvSD.upr = numeric(length = length(dens)))
Pred$dens <- dens
for (i in 1:length(dens)){
# median probability
Pred$FPvSM.mean[i] <- mean(mm3.1$Sol[,1] + dens[i] * mm3.1$Sol[,3])
Pred$FPvSD.mean[i] <- mean(mm3.2$Sol[,2] + dens[i] * mm3.2$Sol[,4])
# lower HPD interval
Pred$FPvSM.lwr[i] <- HPDinterval(mcmc(mm3.1$Sol[,1] + dens[i] * mm3.1$Sol[,3]))[1]
Pred$FPvSD.lwr[i] <- HPDinterval(mcmc(mm3.2$Sol[,2] + dens[i] * mm3.2$Sol[,4]))[1]
# upper HPD interval
Pred$FPvSM.upr[i] <- HPDinterval(mcmc(mm3.1$Sol[,1] + dens[i] * mm3.1$Sol[,3]))[2]
Pred$FPvSD.upr[i] <- HPDinterval(mcmc(mm3.2$Sol[,2] + dens[i] * mm3.2$Sol[,4]))[2]
}Plot log odds
# melt for plotting
Pred$sample <- rownames(Pred)
Pdat <- data.frame(dens = rep(dens,2))
Pdat$FemState <- rep(c("FPvsSM","FPvsSD"), each = length(dens))
Pdat$mean <- reshape::melt(Pred[,c(8,1:2)], id.vars = "sample")[,3]
Pdat$lwr <- reshape::melt(Pred[,c(8,3:4)], id.vars = "sample")[,3]
Pdat$upr <- reshape::melt(Pred[,c(8,5:6)], id.vars = "sample")[,3]
#head(Pdat)
pal <- wes_palette("Zissou1", 12, type = "continuous")
pm3 <- ggplot(data = Pdat, aes(x=dens, y =mean, color = FemState ))+
geom_line(linewidth=1, show.legend = NA)+
theme_bw(base_size = 12)+
geom_ribbon(aes(ymax = upr, ymin = lwr, x=dens, fill=FemState),color=NA, alpha = 0.1)+
ylab(label = "Log-odds ratio of FP vs alternative")+
xlab(label="Log population density")+ #xlim(0,80)+
scale_color_manual(breaks = c("FPvsSM", "FPvsSD"),
values = c(pal[c(5,9)]),
labels = c("SM", "SD"),
name = "Alternative state") +
scale_fill_manual(breaks = c("FPvsSM", "FPvsSD"),
values = c(pal[c(5,9)]),
labels = c("SM", "SD"),
name = "Alternative state") +
geom_hline(yintercept = 0, colour = "gray50", lty = 2) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
pm3
Hypothesis test
n <- nrow(mm3.1$Sol)
Model3.summ <- data.frame(Log_odds = c("FP_vs_SM", "FP_vs_SD"),
PM = c(mean(mm3.1$Sol[,3]), mean(mm3.2$Sol[,4])),
Lwr = c(HPDinterval(mm3.1$Sol[,3])[1],
HPDinterval(mm3.2$Sol[,4])[1]),
Upr = c(HPDinterval(mm3.1$Sol[,3])[2],
HPDinterval(mm3.2$Sol[,4])[2]),
PMCMC = c(sum(mm3.1$Sol[,3] < 0) / n,
sum(mm3.2$Sol[,4] < 0) / n))
Model3.summ %>%
kbl(booktabs =T, align= c(rep('l',1),rep('r', 4))) %>%
kable_styling(latex_options ="striped")| Log_odds | PM | Lwr | Upr | PMCMC |
|---|---|---|---|---|
| FP_vs_SM | 1.601878 | 0.8874501 | 2.281255 | 0.0000000 |
| FP_vs_SD | 1.106785 | 0.3066866 | 1.838425 | 0.0006667 |
OU model
As a second approach to validate these results I modeled the evolution of population density as continuous trait with lineage specific optima. The goal is to assess if polymorphic species have a higher density optimum. First, I need to prepare the data.
dat<-read.csv("../data/processed/glmm/Field_comp_dat.csv", sep = ",", header = T, na.strings = c("", "NA"), colClasses = c(rep("factor", 8), rep("numeric",3)))
# total number or individuals
dat$Total <- dat$Female + dat$Male
# density
dat$dens <- dat$Total / dat$Time * 60
dat_sp <- dat %>% dplyr::select(Taxon,dens,FemState) %>% group_by(Taxon, FemState) %>% summarise(mean = mean(dens))
# re-code colour states as 0 = monomrphic and 1 = polymorphic
dens <- log(dat_sp$mean)
names(dens) <- dat_sp$Taxon
# write female-colour data for correlation with latitude analysis
#write.nexus.data(dens, file = "../data/processed/OU/density.nex", interleaved = T, format = "continuous")
# check taxon labels match
matching <- name.check(tre, dat, dat$Taxon)
# prune tre to match data
denstre <- drop.tip(tre, matching$tree_not_data, rooted = TRUE)
# fix rounding errors to make ultrametric
denstre <- nnls.tree(cophenetic(denstre), denstre, rooted = TRUE, method = "ultrametric")
# write tree for density data
#write.tree(denstre, "../data/processed/OU/density.tree")check_dens <- checkConvergence(list_files = c("~/Dropbox/Doctorado/COEN_Comparative/FP_Comparative/output/OU/Density_OU_run_1.log", "~/Dropbox/Doctorado/COEN_Comparative/FP_Comparative/output/OU/Density_OU_run_2.log"))## Reading in log file 1
## Reading in log file 1
## [1] "Calculating burn-in"
## [1] "Analyzing continuous parameters"## ACHIEVED CONVERGENCE
##
## BURN-IN SET AT 0
##
##
## LOWEST CONTINUOUS PARAMETER ESS
## RUN 1 -> branch_thetas.28. 2562.59
## RUN 2 -> branch_deltas.143. 3078.24
## ## means ESS
## alpha 0.0496481381 7869.217
## branch_deltas.1. -0.0268435117 7753.528
## branch_deltas.2. 0.0016487161 8002.000
## branch_deltas.3. -0.0070707286 7977.832
## branch_deltas.4. -0.0156578315 7315.037
## branch_deltas.5. -0.0052007858 8002.000
## branch_deltas.6. -0.0032129009 7681.852
## branch_deltas.7. 0.0021102836 7493.995
## branch_deltas.8. 0.0040324025 8002.000
## branch_deltas.9. -0.0394167385 7914.330
## branch_deltas.10. 0.0093494481 7921.177
## branch_deltas.11. 0.0269752703 7768.223
## branch_deltas.12. -0.0513211851 8002.000
## branch_deltas.13. 0.0445869251 7874.757
## branch_deltas.14. 0.0043364794 7644.493
## branch_deltas.15. -0.0666836463 8002.000
## branch_deltas.16. -0.0146889002 7555.851
## branch_deltas.17. -0.0241626231 7539.747
## branch_deltas.18. 0.0158157950 7520.988
## branch_deltas.19. -0.0001291965 8002.000
## branch_deltas.20. 0.0052414194 7583.246
## branch_deltas.21. -0.0316391272 8002.000
## branch_deltas.22. -0.0200172731 8002.000
## branch_deltas.23. -0.0013132084 7791.797
## branch_deltas.24. 0.0135167000 7432.235
## branch_deltas.25. -0.0007979157 7047.455
## branch_deltas.26. 0.0088906029 8002.000
## branch_deltas.27. 0.0078478107 7587.088
## branch_deltas.28. -0.0169256250 8002.000
## branch_deltas.29. 0.0007229185 8002.000
## branch_deltas.30. 0.0049142480 7676.382
## branch_deltas.31. 0.0170615927 8002.000
## branch_deltas.32. -0.0031107556 8002.000
## branch_deltas.33. -0.0140581240 7209.652
## branch_deltas.34. 0.0084738709 6820.862
## branch_deltas.35. -0.0144992743 7829.821
## branch_deltas.36. 0.0129485098 8002.000
## branch_deltas.37. 0.0118790119 8002.000
## branch_deltas.38. -0.0038866362 7556.995
## branch_deltas.39. 0.0003138962 7875.404
## branch_deltas.40. 0.0057697133 7828.897
## branch_deltas.41. 0.0011351100 8002.000
## branch_deltas.42. -0.0036577220 7770.444
## branch_deltas.43. 0.0439435623 8002.000
## branch_deltas.44. 0.0360615286 7781.688
## branch_deltas.45. 0.0015621406 7751.823
## branch_deltas.46. -0.0057836532 7918.374
## branch_deltas.47. -0.0208275693 7785.099
## branch_deltas.48. -0.0503175931 7975.536
## branch_deltas.49. 0.0412045490 8002.000
## branch_deltas.50. -0.0214639847 7898.977
## branch_deltas.51. -0.0030940032 7701.216
## branch_deltas.52. -0.0193172183 7335.866
## branch_deltas.53. 0.0152327803 7839.749
## branch_deltas.54. 0.0151257170 7289.534
## branch_deltas.55. 0.0113807922 7115.443
## branch_deltas.56. 0.0205557770 7789.612
## branch_deltas.57. -0.0062092512 8002.000
## branch_deltas.58. -0.0270001680 7388.448
## branch_deltas.59. -0.0002957047 7801.793
## branch_deltas.60. 0.0351200742 8002.000
## branch_deltas.61. -0.0006721918 7470.609
## branch_deltas.62. 0.0156907011 7897.920
## branch_deltas.63. 0.0343732170 7834.689
## branch_deltas.64. -0.0282755599 7773.689
## branch_deltas.65. 0.0034772023 7544.698
## branch_deltas.66. 0.0116727212 7979.055
## branch_deltas.67. 0.0130975388 7656.321
## branch_deltas.68. -0.0262085499 7974.367
## branch_deltas.69. -0.0021937100 7528.216
## branch_deltas.70. 0.0283925858 7767.820
## branch_deltas.71. -0.0215767174 7765.189
## branch_deltas.72. 0.0122136555 8002.000
## branch_deltas.73. 0.0155291082 7359.596
## branch_deltas.74. -0.0143250453 8002.000
## branch_deltas.75. 0.0120361979 7900.090
## branch_deltas.76. 0.0053136423 7979.036
## branch_deltas.77. 0.0104257060 7513.907
## branch_deltas.78. -0.0646807313 7439.884
## branch_deltas.79. 0.0256336704 7920.802
## branch_deltas.80. -0.0052428056 7352.712
## branch_deltas.81. -0.0038594161 7502.789
## branch_deltas.82. -0.0014260041 8002.000
## branch_deltas.83. -0.0101404149 7741.612
## branch_deltas.84. -0.0131413592 7835.328
## branch_deltas.85. -0.0224430787 8002.000
## branch_deltas.86. 0.0239048034 7869.459
## branch_deltas.87. -0.0281213460 7980.986
## branch_deltas.88. -0.0544152380 7880.354
## branch_deltas.89. 0.0182723669 7817.706
## branch_deltas.90. 0.0588713046 7697.451
## branch_deltas.91. -0.0038766696 7982.319
## branch_deltas.92. 0.0010152117 8002.000
## branch_deltas.93. 0.0080282806 8002.000
## branch_deltas.94. -0.0001223668 8002.000
## branch_deltas.95. -0.0088961652 7856.114
## branch_deltas.96. -0.0339586465 7888.649
## branch_deltas.97. -0.0089307730 8002.000
## branch_deltas.98. 0.0059358622 7226.906
## branch_deltas.99. 0.0285876390 8002.000
## branch_deltas.100. 0.0022017605 7771.011
## branch_deltas.101. 0.0215680559 7869.456
## branch_deltas.102. 0.0297479531 7390.099
## branch_deltas.103. 0.0298646043 7683.836
## branch_deltas.104. 0.0236298792 7001.679
## branch_deltas.105. 0.0601286015 7742.875
## branch_deltas.106. 0.0697204399 7986.387
## branch_deltas.107. -0.0547489865 7337.843
## branch_deltas.108. 0.0147999202 7633.498
## branch_deltas.109. 0.0018227280 7982.618
## branch_deltas.110. 0.0489369598 8002.000
## branch_deltas.111. 0.0882123332 8002.000
## branch_deltas.112. 0.0402146541 8002.000
## branch_deltas.113. -0.0004634346 7570.295
## branch_deltas.114. -0.0269898603 7816.808
## branch_deltas.115. -0.0278199378 7749.273
## branch_deltas.116. -0.0012087482 7258.827
## branch_deltas.117. -0.0230850640 8002.000
## branch_deltas.118. -0.0236376538 7126.433
## branch_deltas.119. -0.0555233890 7738.951
## branch_deltas.120. -0.0507869160 7678.084
## branch_deltas.121. -0.0164191986 7721.119
## branch_deltas.122. 0.0001609194 7840.813
## branch_deltas.123. 0.0292983879 7902.524
## branch_deltas.124. 0.1078269968 7825.777
## branch_deltas.125. 0.0234030956 7495.439
## branch_deltas.126. 0.1015557144 7104.372
## branch_deltas.127. -0.0038329505 7858.789
## branch_deltas.128. 0.0273533628 7901.303
## branch_deltas.129. -0.0139813110 7978.543
## branch_deltas.130. -0.0321066403 7542.642
## branch_deltas.131. -0.0039415643 7954.727
## branch_deltas.132. -0.0025038555 8002.000
## branch_deltas.133. -0.0037091949 7457.270
## branch_deltas.134. 0.0077130577 8002.000
## branch_deltas.135. 0.0048043651 7882.438
## branch_deltas.136. -0.0139542757 7779.046
## branch_deltas.137. 0.0121950924 8002.000
## branch_deltas.138. 0.0013259569 8002.000
## branch_deltas.139. -0.0019252046 7839.038
## branch_deltas.140. 0.0045206696 7933.400
## branch_deltas.141. 0.0232650091 8002.000
## branch_deltas.142. 0.0178230935 7621.728
## branch_deltas.143. -0.0521030222 7079.241
## branch_deltas.144. 0.0185923421 7511.042
## branch_deltas.145. 0.0397972909 7924.385
## branch_deltas.146. 0.0096853984 7889.302
## branch_deltas.147. -0.0022835391 7248.933
## branch_deltas.148. -0.0363424813 7595.900
## branch_deltas.149. -0.0481719306 7984.136
## branch_deltas.150. -0.1038672651 8002.000
## branch_deltas.151. -0.0042613728 7513.321
## branch_deltas.152. 0.0390384225 6631.684
## branch_deltas.153. 0.0255392901 8002.000
## branch_deltas.154. -0.0349454405 7653.139
## branch_deltas.155. -0.0097135325 7888.732
## branch_deltas.156. -0.0290438648 7442.619
## branch_deltas.157. -0.0032134456 7358.682
## branch_deltas.158. -0.0194978982 7958.094
## branch_deltas.159. -0.0252612383 8002.000
## branch_deltas.160. -0.0227791917 7929.869
## branch_deltas.161. -0.0169157702 7756.914
## branch_deltas.162. -0.0479110273 8002.000
## branch_deltas.163. -0.0816230718 7709.739
## branch_deltas.164. -0.1063414919 7779.802
## branch_theta_shift.1. 0.1178455386 7458.586
## branch_theta_shift.2. 0.1128467883 7559.645
## branch_theta_shift.3. 0.1132216946 7754.475
## branch_theta_shift.4. 0.1120969758 7818.170
## branch_theta_shift.5. 0.1149712572 8002.000
## branch_theta_shift.6. 0.1218445389 7686.469
## branch_theta_shift.7. 0.1138465384 7566.463
## branch_theta_shift.8. 0.1233441640 7643.836
## branch_theta_shift.9. 0.1320919770 8002.000
## branch_theta_shift.10. 0.1209697576 7613.988
## branch_theta_shift.11. 0.1270932267 6306.678
## branch_theta_shift.12. 0.1413396651 7807.061
## branch_theta_shift.13. 0.1269682579 7770.788
## branch_theta_shift.14. 0.1099725069 7648.604
## branch_theta_shift.15. 0.1529617596 7703.232
## branch_theta_shift.16. 0.1070982254 7224.661
## branch_theta_shift.17. 0.1213446638 7659.180
## branch_theta_shift.18. 0.1160959760 7821.453
## branch_theta_shift.19. 0.1195951012 7773.902
## branch_theta_shift.20. 0.1225943514 8002.000
## branch_theta_shift.21. 0.1184703824 7560.157
## branch_theta_shift.22. 0.1118470382 7612.961
## branch_theta_shift.23. 0.1233441640 7260.702
## branch_theta_shift.24. 0.1279680080 7861.087
## branch_theta_shift.25. 0.1079730067 8002.000
## branch_theta_shift.26. 0.1125968508 7782.980
## branch_theta_shift.27. 0.1137215696 7974.601
## branch_theta_shift.28. 0.1147213197 7079.080
## branch_theta_shift.29. 0.1148462884 8002.000
## branch_theta_shift.30. 0.1077230692 7828.252
## branch_theta_shift.31. 0.1188452887 7809.276
## branch_theta_shift.32. 0.1134716321 7994.082
## branch_theta_shift.33. 0.1205948513 6989.637
## branch_theta_shift.34. 0.1203449138 7484.479
## branch_theta_shift.35. 0.1164708823 7775.601
## branch_theta_shift.36. 0.1253436641 7859.532
## branch_theta_shift.37. 0.1188452887 7566.169
## branch_theta_shift.38. 0.1199700075 7566.886
## branch_theta_shift.39. 0.1155961010 7985.022
## branch_theta_shift.40. 0.1245938515 8002.000
## branch_theta_shift.41. 0.1220944764 8002.000
## branch_theta_shift.42. 0.1160959760 7626.658
## branch_theta_shift.43. 0.1337165709 7917.726
## branch_theta_shift.44. 0.1219695076 7877.007
## branch_theta_shift.45. 0.1108472882 8002.000
## branch_theta_shift.46. 0.1217195701 7921.861
## branch_theta_shift.47. 0.1233441640 8002.000
## branch_theta_shift.48. 0.1390902274 7918.390
## branch_theta_shift.49. 0.1360909773 8002.000
## branch_theta_shift.50. 0.1195951012 7804.310
## branch_theta_shift.51. 0.1142214446 7400.315
## branch_theta_shift.52. 0.1168457886 7751.207
## branch_theta_shift.53. 0.1132216946 7289.908
## branch_theta_shift.54. 0.1185953512 7765.101
## branch_theta_shift.55. 0.1192201950 7938.931
## branch_theta_shift.56. 0.1228442889 7639.992
## branch_theta_shift.57. 0.1148462884 7867.897
## branch_theta_shift.58. 0.1243439140 7678.633
## branch_theta_shift.59. 0.1084728818 7792.273
## branch_theta_shift.60. 0.1319670082 7322.350
## branch_theta_shift.61. 0.1172206948 7969.845
## branch_theta_shift.62. 0.1180954761 7771.543
## branch_theta_shift.63. 0.1268432892 7144.842
## branch_theta_shift.64. 0.1234691327 7943.871
## branch_theta_shift.65. 0.1258435391 7757.505
## branch_theta_shift.66. 0.1295926018 7544.093
## branch_theta_shift.67. 0.1247188203 7850.404
## branch_theta_shift.68. 0.1324668833 7887.085
## branch_theta_shift.69. 0.1289677581 7915.643
## branch_theta_shift.70. 0.1283429143 7748.265
## branch_theta_shift.71. 0.1189702574 7929.880
## branch_theta_shift.72. 0.1162209448 8002.000
## branch_theta_shift.73. 0.1223444139 7396.717
## branch_theta_shift.74. 0.1140964759 8002.000
## branch_theta_shift.75. 0.1160959760 7988.690
## branch_theta_shift.76. 0.1239690077 7371.477
## branch_theta_shift.77. 0.1265933517 7799.926
## branch_theta_shift.78. 0.1495876031 7501.605
## branch_theta_shift.79. 0.1232191952 8002.000
## branch_theta_shift.80. 0.1179705074 8002.000
## branch_theta_shift.81. 0.1197200700 7613.580
## branch_theta_shift.82. 0.1169707573 7799.942
## branch_theta_shift.83. 0.1205948513 7942.513
## branch_theta_shift.84. 0.1133466633 7588.765
## branch_theta_shift.85. 0.1153461635 8002.000
## branch_theta_shift.86. 0.1203449138 7701.591
## branch_theta_shift.87. 0.1180954761 8002.000
## branch_theta_shift.88. 0.1404648838 7994.079
## branch_theta_shift.89. 0.1214696326 7678.968
## branch_theta_shift.90. 0.1482129468 7791.585
## branch_theta_shift.91. 0.0934766308 8002.000
## branch_theta_shift.92. 0.1062234441 7586.647
## branch_theta_shift.93. 0.1080979755 7691.403
## branch_theta_shift.94. 0.0884778805 7470.508
## branch_theta_shift.95. 0.0903524119 8002.000
## branch_theta_shift.96. 0.1342164459 8002.000
## branch_theta_shift.97. 0.1177205699 7737.333
## branch_theta_shift.98. 0.1122219445 7790.348
## branch_theta_shift.99. 0.1189702574 7901.735
## branch_theta_shift.100. 0.1114721320 7609.096
## branch_theta_shift.101. 0.1089727568 7137.670
## branch_theta_shift.102. 0.1148462884 7544.755
## branch_theta_shift.103. 0.1215946013 7325.757
## branch_theta_shift.104. 0.1178455386 7226.833
## branch_theta_shift.105. 0.1447138215 7949.559
## branch_theta_shift.106. 0.1535866033 7714.708
## branch_theta_shift.107. 0.1424643839 7721.778
## branch_theta_shift.108. 0.0954761310 7686.808
## branch_theta_shift.109. 0.0891027243 7905.425
## branch_theta_shift.110. 0.1332166958 8002.000
## branch_theta_shift.111. 0.1719570107 7784.502
## branch_theta_shift.112. 0.1200949763 8002.000
## branch_theta_shift.113. 0.1069732567 7708.741
## branch_theta_shift.114. 0.1189702574 7458.680
## branch_theta_shift.115. 0.1152211947 7918.254
## branch_theta_shift.116. 0.1158460385 8002.000
## branch_theta_shift.117. 0.1140964759 8002.000
## branch_theta_shift.118. 0.1110972257 7421.754
## branch_theta_shift.119. 0.1372156961 7820.348
## branch_theta_shift.120. 0.1265933517 7533.571
## branch_theta_shift.121. 0.1029742564 8002.000
## branch_theta_shift.122. 0.1244688828 7365.581
## branch_theta_shift.123. 0.1248437891 7921.129
## branch_theta_shift.124. 0.1905773557 7619.864
## branch_theta_shift.125. 0.1028492877 7832.209
## branch_theta_shift.126. 0.2009497626 7612.429
## branch_theta_shift.127. 0.1082229443 8002.000
## branch_theta_shift.128. 0.1243439140 7970.311
## branch_theta_shift.129. 0.1159710072 8002.000
## branch_theta_shift.130. 0.1229692577 8002.000
## branch_theta_shift.131. 0.0998500375 7572.178
## branch_theta_shift.132. 0.0974756311 8002.000
## branch_theta_shift.133. 0.1002249438 7653.361
## branch_theta_shift.134. 0.0974756311 8002.000
## branch_theta_shift.135. 0.1089727568 7861.968
## branch_theta_shift.136. 0.1102224444 8002.000
## branch_theta_shift.137. 0.1063484129 8002.000
## branch_theta_shift.138. 0.0956010997 8002.000
## branch_theta_shift.139. 0.0986003499 7443.908
## branch_theta_shift.140. 0.0834791302 7179.300
## branch_theta_shift.141. 0.1239690077 7525.031
## branch_theta_shift.142. 0.0927268183 7435.648
## branch_theta_shift.143. 0.1344663834 7821.005
## branch_theta_shift.144. 0.1252186953 7584.307
## branch_theta_shift.145. 0.1317170707 7764.028
## branch_theta_shift.146. 0.1049737566 7970.202
## branch_theta_shift.147. 0.1024743814 7879.603
## branch_theta_shift.148. 0.1168457886 7972.542
## branch_theta_shift.149. 0.1287178205 8002.000
## branch_theta_shift.150. 0.1923269183 7988.727
## branch_theta_shift.151. 0.1087228193 7748.813
## branch_theta_shift.152. 0.1312171957 7770.426
## branch_theta_shift.153. 0.1079730067 7322.166
## branch_theta_shift.154. 0.1174706323 8002.000
## branch_theta_shift.155. 0.0901024744 7995.418
## branch_theta_shift.156. 0.1210947263 7719.204
## branch_theta_shift.157. 0.1155961010 7523.249
## branch_theta_shift.158. 0.1102224444 7435.826
## branch_theta_shift.159. 0.1139715071 7797.242
## branch_theta_shift.160. 0.1178455386 7868.305
## branch_theta_shift.161. 0.1120969758 8002.000
## branch_theta_shift.162. 0.1342164459 7862.071
## branch_theta_shift.163. 0.1755811047 7339.127
## branch_theta_shift.164. 0.2121969508 7461.095
## branch_thetas.1. 1.6599936342 7941.324
## branch_thetas.2. 1.6236590824 6614.710
## branch_thetas.3. 1.5921604238 6588.728
## branch_thetas.4. 1.5835733214 7873.731
## branch_thetas.5. 1.6084641108 7917.812
## branch_thetas.6. 1.5877406301 7898.046
## branch_thetas.7. 1.5930638188 8002.000
## branch_thetas.8. 1.5691555394 7956.780
## branch_thetas.9. 1.5257063828 8002.000
## branch_thetas.10. 1.7977283957 7870.982
## branch_thetas.11. 1.8153542178 8002.000
## branch_thetas.12. 1.6937579707 7769.261
## branch_thetas.13. 1.7896660928 7343.216
## branch_thetas.14. 1.5760984865 7687.388
## branch_thetas.15. 1.5050783959 7556.812
## branch_thetas.16. 1.5666190184 7453.870
## branch_thetas.17. 1.5668306711 7574.819
## branch_thetas.18. 1.6466063836 7550.298
## branch_thetas.19. 1.6492537532 7786.260
## branch_thetas.20. 1.6546243666 7633.942
## branch_thetas.21. 1.4998492934 7803.510
## branch_thetas.22. 1.5114711458 7670.781
## branch_thetas.23. 1.7985215834 7753.330
## branch_thetas.24. 1.8133515042 7245.711
## branch_thetas.25. 1.7783673025 6713.198
## branch_thetas.26. 1.8015768810 7841.560
## branch_thetas.27. 1.8005340810 7354.279
## branch_thetas.28. 1.7496113134 6008.617
## branch_thetas.29. 1.7720642250 7557.546
## branch_thetas.30. 1.7762555424 6690.544
## branch_thetas.31. 1.8058650987 7648.743
## branch_thetas.32. 1.7819835873 7308.372
## branch_thetas.33. 1.7185027927 7659.498
## branch_thetas.34. 1.7410348286 8002.000
## branch_thetas.35. 1.7320429925 8002.000
## branch_thetas.36. 1.8189507538 7909.489
## branch_thetas.37. 1.8178812720 8002.000
## branch_thetas.38. 1.7748708998 7561.501
## branch_thetas.39. 1.7790713894 7497.414
## branch_thetas.40. 2.1426555567 7995.933
## branch_thetas.41. 2.1381818717 8002.000
## branch_thetas.42. 2.1333890630 8002.000
## branch_thetas.43. 2.1515310517 8002.000
## branch_thetas.44. 2.0194027639 7776.395
## branch_thetas.45. 1.9341164954 7908.844
## branch_thetas.46. 1.8476096328 7842.611
## branch_thetas.47. 1.8094806833 7946.754
## branch_thetas.48. 1.7787819035 8002.000
## branch_thetas.49. 1.8703040097 7896.327
## branch_thetas.50. 1.8007572401 7493.396
## branch_thetas.51. 1.8191271682 7448.487
## branch_thetas.52. 1.8294303959 7904.325
## branch_thetas.53. 1.8639803963 7876.370
## branch_thetas.54. 2.1688470477 7865.993
## branch_thetas.55. 2.1651021456 7972.708
## branch_thetas.56. 2.1253401810 7905.991
## branch_thetas.57. 1.9722363896 7926.518
## branch_thetas.58. 1.9514454695 7821.814
## branch_thetas.59. 2.1863778834 7876.496
## branch_thetas.60. 2.2516582476 7781.879
## branch_thetas.61. 2.2158659995 7720.140
## branch_thetas.62. 2.2084823534 7814.917
## branch_thetas.63. 2.2487329288 7772.695
## branch_thetas.64. 2.1882858908 7874.934
## branch_thetas.65. 2.2486263385 7014.369
## branch_thetas.66. 2.2568218791 7504.802
## branch_thetas.67. 2.1219485059 7943.787
## branch_thetas.68. 2.0397530345 7934.754
## branch_thetas.69. 2.0637678622 7932.900
## branch_thetas.70. 2.1283127770 7987.878
## branch_thetas.71. 1.9879218712 7835.995
## branch_thetas.72. 2.0296181374 7839.705
## branch_thetas.73. 2.0339487793 8002.000
## branch_thetas.74. 2.0040946150 7844.046
## branch_thetas.75. 2.0764069984 7809.647
## branch_thetas.76. 2.0879568171 7769.319
## branch_thetas.77. 2.0930688565 7906.650
## branch_thetas.78. 1.8582821987 6960.080
## branch_thetas.79. 1.9725014334 7566.071
## branch_thetas.80. 1.9416249642 7423.182
## branch_thetas.81. 1.9116404496 7701.159
## branch_thetas.82. 1.9140738620 7723.704
## branch_thetas.83. 1.9185007946 7667.743
## branch_thetas.84. 1.9154998422 7718.396
## branch_thetas.85. 1.9286412242 7786.345
## branch_thetas.86. 1.9468677655 7651.945
## branch_thetas.87. 1.9229629391 7724.982
## branch_thetas.88. 1.9510842968 7625.557
## branch_thetas.89. 2.0826431653 7753.236
## branch_thetas.90. 2.0643708123 7741.759
## branch_thetas.91. 2.0054995320 7643.518
## branch_thetas.92. 2.0184196680 8002.000
## branch_thetas.93. 2.0174044510 8002.000
## branch_thetas.94. 2.0093761676 7838.735
## branch_thetas.95. 2.0094985707 7867.315
## branch_thetas.96. 2.0659615598 8002.000
## branch_thetas.97. 2.0999202001 7989.600
## branch_thetas.98. 2.1088509681 8002.000
## branch_thetas.99. 2.2451491231 7857.389
## branch_thetas.100. 2.2165614859 7879.821
## branch_thetas.101. 2.2143597272 7860.950
## branch_thetas.102. 2.1927916562 7868.885
## branch_thetas.103. 2.2165381767 7705.466
## branch_thetas.104. 2.1866735881 7952.570
## branch_thetas.105. 2.1630436905 7981.219
## branch_thetas.106. 2.1029150955 8002.000
## branch_thetas.107. 1.9784456463 7923.683
## branch_thetas.108. 2.0331946563 8002.000
## branch_thetas.109. 2.0183947497 7839.480
## branch_thetas.110. 2.1537213430 7949.797
## branch_thetas.111. 2.1047843869 7928.773
## branch_thetas.112. 2.0165720145 7917.977
## branch_thetas.113. 1.8487476266 7931.204
## branch_thetas.114. 1.8222212085 7668.143
## branch_thetas.115. 1.8492110616 7915.078
## branch_thetas.116. 1.8290994764 7985.603
## branch_thetas.117. 1.8303082247 7837.047
## branch_thetas.118. 1.8533933160 7784.142
## branch_thetas.119. 1.8770309580 7934.953
## branch_thetas.120. 1.9325543608 7998.404
## branch_thetas.121. 1.9833412621 7773.792
## branch_thetas.122. 2.1370467489 8002.000
## branch_thetas.123. 2.1368858433 7983.372
## branch_thetas.124. 2.1075874840 7934.975
## branch_thetas.125. 1.9997604604 7892.778
## branch_thetas.126. 1.9763573657 7990.487
## branch_thetas.127. 1.7787575068 7568.680
## branch_thetas.128. 1.8060022626 7950.350
## branch_thetas.129. 1.7325609387 7790.825
## branch_thetas.130. 1.7465422697 7808.373
## branch_thetas.131. 1.7786489123 7770.256
## branch_thetas.132. 1.7825904624 7709.949
## branch_thetas.133. 1.7850943280 7318.617
## branch_thetas.134. 1.7888035110 7269.530
## branch_thetas.135. 1.7713413191 7638.218
## branch_thetas.136. 1.7665369347 7388.609
## branch_thetas.137. 1.7926862778 7487.690
## branch_thetas.138. 1.7804912112 6956.947
## branch_thetas.139. 1.7791652407 7371.630
## branch_thetas.140. 1.7810904483 7453.889
## branch_thetas.141. 1.7998347978 7226.072
## branch_thetas.142. 1.7765697904 7645.109
## branch_thetas.143. 1.5314884230 7410.037
## branch_thetas.144. 1.6493829511 7785.248
## branch_thetas.145. 1.6307906026 7817.881
## branch_thetas.146. 1.5909933091 7608.136
## branch_thetas.147. 1.5813079112 7838.172
## branch_thetas.148. 1.5835914566 8002.000
## branch_thetas.149. 1.5717620310 7815.363
## branch_thetas.150. 1.6199339610 8001.436
## branch_thetas.151. 1.7450791475 7476.118
## branch_thetas.152. 1.7883789477 7558.364
## branch_thetas.153. 1.7493405072 7485.182
## branch_thetas.154. 1.7238012294 7178.063
## branch_thetas.155. 1.7587466838 7137.941
## branch_thetas.156. 1.5651231369 7908.403
## branch_thetas.157. 1.5909535583 8001.279
## branch_thetas.158. 1.5941670058 7880.407
## branch_thetas.159. 1.6136648951 7858.241
## branch_thetas.160. 1.5992311555 7062.937
## branch_thetas.161. 1.6220103542 6828.543
## branch_thetas.162. 1.6389261335 7928.294
## branch_thetas.163. 1.6868371493 7863.340
## branch_thetas.164. 1.7684602051 7099.159
## num_theta_changes 19.8030492377 7106.502
## sigma2 0.1553742747 7898.596
## theta_root 1.8748016593 7334.919
# read the annotated tree
tree <- read.beast("../output/OU/Density_OU_MAP.tre")
#tree <- read.beast("../output/OU/Density_OU_10_shifts_MAP.tre")
#tree <- read.beast("../output/OU/Density_OU_40_shifts_MAP.tre")
d <- data.frame(dat_sp$FemState)
d$dat_sp.FemState <- as.factor(d$dat_sp.FemState)
levels(d$dat_sp.FemState) <- c("0","1", "0")
rownames(d) <- gsub("_", " ", dat_sp$Taxon)
tip_labels <- gsub("_", " ", tree@phylo$tip.label)
tree@phylo$tip.label <- tip_labels
pal <- wes_palette("Zissou1", 12, type = "continuous")
p1 <- ggtree(tree, aes(color=branch_thetas), size = 0.5) +
geom_tiplab(size =1.5, colour = "black", offset = 10) +
scale_colour_gradientn(colours = pal, name = "Log(density)") +
xlim(0,210)
g1 <- gheatmap(p = p1, data = d, offset=0.5, width=0.05,
colnames=FALSE, legend_title="Colour state", color = "black") +
scale_fill_manual(values = c("white", "grey10"),
name = "Female colour \nstate",
labels = c("FM", "FP"))
g1
#ggsave(g1, filename = "../../Writing/Figures//Density_OU_branch_thetas.pdf", device = "pdf", dpi = 600, units = "mm", width = 90, height = 115)
#ggsave(g1, filename = "../../Writing/Figures//Density_OU_40_branch_thetas.pdf", device = "pdf", dpi = 600, units = "mm", width = 90, height = 115)
#ggsave(g1, filename = "../../Writing/Figures//Density_OU_0_branch_thetas.pdf", device = "pdf", dpi = 600, units = "mm", width = 90, height = 115)Effects of OSR on the occurrence of FPs
Run the mode
inv.phylo<-inverseA(ecotre2,nodes="TIPS",scale=TRUE)
#residual covariance matrix 1/J*(I+J), I and J are identity and unit matrices respectively
IJ <- (1/3) * (diag(2) + matrix(1, 2, 2))
# kronecker prior for fixed effect close to being flat for the 2 way:
#(0 vs.1 , 1 vs. 2 and 0 vs. 2) marginal probabilities within each fixed effect
prDem2 = list(
B = list(mu = rep(0, 4), V = diag(4) * (1.7 + pi ^ 2 / 3)),
R = list(V = IJ, fix = 1),
G = list(
G1 = list(
V = 1,
nu = 1000,
alpha.mu = rep(0, 1),
alpha.V = diag(1)
),
G2 = list(
V = 1,
nu = 1000,
alpha.mu = rep(0, 1),
alpha.V = diag(1)
)
)
)
#The model for the effect of OSR
mm4.1 <- MCMCglmm(
relevel(FemState, ref = "0") ~ trait + log(OSR):trait,
random = ~ Site + Taxon,
rcov = ~ us(trait):units,
ginverse = list(Taxon = inv.phylo$Ainv),
family = "categorical",
data = ecodat2,
prior = prDem2,
pl = TRUE,
slice = TRUE,
nitt = 330000,
burnin = 30000,
thin = 100,
verbose = F
)
mm4.2 <- MCMCglmm(
relevel(FemState, ref = "2") ~ trait + log(OSR):trait,
random = ~ Site + Taxon,
rcov = ~ us(trait):units,
ginverse = list(Taxon = inv.phylo$Ainv),
family = "categorical",
data = ecodat2,
prior = prDem2,
pl = TRUE,
slice = TRUE,
nitt = 330000,
burnin = 30000,
thin = 100,
verbose = F
)Get estimates in response scale
# rescale estimates assuming residual covariance matrix of zero
Delta <- cbind(c(-1, 1, 0), c(-1, 0, 1))
c2 <- (16 * sqrt(3)/(15 * pi))^2
D <- ginv(Delta %*% t(Delta)) %*% Delta
OSR <- log(ecodat2$OSR)
# make the estimate data frame
Pred <- data.frame(SM.mean = numeric(length = length(OSR)),
FP.mean = numeric(length = length(OSR)),
SD.mean = numeric(length = length(OSR)),
SM.lwr = numeric(length = length(OSR)),
FP.lwr = numeric(length = length(OSR)),
SD.lwr = numeric(length = length(OSR)),
SM.upr = numeric(length = length(OSR)),
FP.upr = numeric(length = length(OSR)),
SD.upr = numeric(length = length(OSR)))
Pred$OSR <- OSR
Int <- data.frame(SM = numeric(length = 2*nrow(mm4.1$Sol)),
FP = numeric(length = 2*nrow(mm4.1$Sol)),
SD = numeric(length = 2*nrow(mm4.1$Sol)))
for (i in 1:length(OSR)){
# rescale the estimates as if the residual covariance matrix was zero
Int1 <- t(apply((mm4.1$Sol[, 1:2]+OSR[i]*mm4.1$Sol[, 3:4]), 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Int2 <- t(apply((mm4.2$Sol[, 1:2]+OSR[i]*mm4.2$Sol[, 3:4]), 1, function(x) {
D %*% (x/sqrt(1 + c2 * diag(IJ))
)
}))
Int <- rbind(Int1, Int2[,c(2,3,1)])
# median probability
Pred$SM.mean[i] <- mean(mcmc(exp(Int[,1])/rowSums(exp(Int))))
Pred$FP.mean[i] <- mean(mcmc(exp(Int[,2])/rowSums(exp(Int))))
Pred$SD.mean[i] <- mean(mcmc(exp(Int[,3])/rowSums(exp(Int))))
# lower HPD interval
Pred$SM.lwr[i] <- HPDinterval(mcmc(exp(Int[,1])/rowSums(exp(Int))))[1]
Pred$FP.lwr[i] <- HPDinterval(mcmc(exp(Int[,2])/rowSums(exp(Int))))[1]
Pred$SD.lwr[i] <- HPDinterval(mcmc(exp(Int[,3])/rowSums(exp(Int))))[1]
# upper HPD interval
Pred$SM.upr[i] <- HPDinterval(mcmc(exp(Int[,1])/rowSums(exp(Int))))[2]
Pred$FP.upr[i] <- HPDinterval(mcmc(exp(Int[,2])/rowSums(exp(Int))))[2]
Pred$SD.upr[i] <- HPDinterval(mcmc(exp(Int[,3])/rowSums(exp(Int))))[2]
}Plot predicted
# melt for plotting
Pred$sample<-rownames(Pred)
Pdat<-data.frame(OSR = rep(OSR,3))
Pdat$FemState<-rep(c("SM","FP", "SD"), each = length(OSR))
Pdat$mean<-reshape::melt(Pred[,c(11,1:3)], id.vars = "sample")[,3]
Pdat$lwr<-reshape::melt(Pred[,c(11,4:6)], id.vars = "sample")[,3]
Pdat$upr<-reshape::melt(Pred[,c(11,7:9)], id.vars = "sample")[,3]
#head(Pdat)
# get species data points
Mdata<-data.frame(taxon = ecodat2$Taxon,
OSR = log(ecodat2$OSR),
SM = numeric(length(ecodat2$Taxon)),
FP = numeric(length(ecodat2$Taxon)),
SD = numeric(length(ecodat2$Taxon)))
for (i in 1:nrow(Mdata)){
if(ecodat2$FemState[i] == 0){
Mdata$SM[i] <- 1}
else{
if(ecodat2$FemState[i] == 1){
Mdata$FP[i] <-1}
else{
Mdata$SD[i] <-1
}
}
}
#head(Mdata)
# melt for plotting
MdataPlot<-reshape::melt(Mdata, id.vars = c("taxon", "OSR"))
colnames(MdataPlot)[3:4]<-c("FemState", "mean")
pal <- wes_palette("Zissou1", 12, type = "continuous")
p4.2<-ggplot(data = Pdat, aes(x=OSR, y =mean, color = FemState ))+
geom_line(linewidth=0.5, show.legend = NA)+
theme_bw(base_size = 8)+
geom_ribbon(aes(ymax = upr, ymin = lwr, x=OSR, fill=FemState),color=NA, alpha = 0.1)+
geom_point(data = MdataPlot, aes(color = FemState, fill = FemState), alpha = 0.3 , size = 1, position = position_jitter(height = 0,) )+
ylab(label = "Probability of \nsex-related colour state")+
xlab(label="Log adult OSR")+ # xlim(0,80)+
scale_color_manual(breaks = c("FP", "SD", "SM"),
values = c(pal[c(11,7,1)]),
name = "Colour state") +
scale_fill_manual(breaks = c("FP", "SD", "SM"),
values = c(pal[c(11,7,1)]),
name = "Colour state") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p4.2
Get log odds estimates
OSR <- log(ecodat2$OSR)
# make the estimate data frame
Pred <- data.frame(FPvSM.mean = numeric(length = length(OSR)),
FPvSD.mean = numeric(length = length(OSR)),
FPvSM.lwr = numeric(length = length(OSR)),
FPvSD.lwr = numeric(length = length(OSR)),
FPvSM.upr = numeric(length = length(OSR)),
FPvSD.upr = numeric(length = length(OSR)))
Pred$OSR <- OSR
Int<-data.frame(MA = numeric(length = nrow(mm4.1$Sol)),
FP = numeric(length = nrow(mm4.1$Sol)),
MH = numeric(length = nrow(mm4.1$Sol)))
for (i in 1:length(OSR)){
# median probability
Pred$FPvSM.mean[i] <- mean(mm4.1$Sol[,1] + OSR[i] * mm4.1$Sol[,3])
Pred$FPvSD.mean[i] <- mean(mm4.2$Sol[,2] + OSR[i] * mm4.2$Sol[,4])
# lower HPD interval
Pred$FPvSM.lwr[i] <- HPDinterval(mcmc(mm4.1$Sol[,1] + OSR[i] * mm4.1$Sol[,3]))[1]
Pred$FPvSD.lwr[i] <- HPDinterval(mcmc(mm4.2$Sol[,2] + OSR[i] * mm4.2$Sol[,4]))[1]
# upper HPD interval
Pred$FPvSM.upr[i] <- HPDinterval(mcmc(mm4.1$Sol[,1] + OSR[i] * mm4.1$Sol[,3]))[2]
Pred$FPvSD.upr[i] <- HPDinterval(mcmc(mm4.2$Sol[,2] + OSR[i] * mm4.2$Sol[,4]))[2]
}Plot log odds
# melt for plotting
Pred$sample<-rownames(Pred)
Pdat<-data.frame(OSR = rep(OSR,2))
Pdat$FemState<-rep(c("FPvsSM","FPvsSD"), each = length(OSR))
Pdat$mean<-reshape::melt(Pred[,c(8,1:2)], id.vars = "sample")[,3]
Pdat$lwr<-reshape::melt(Pred[,c(8,3:4)], id.vars = "sample")[,3]
Pdat$upr<-reshape::melt(Pred[,c(8,5:6)], id.vars = "sample")[,3]
#head(Pdat)
pal <- wes_palette("Zissou1", 12, type = "continuous")
pm4 <-ggplot(data = Pdat, aes(x=OSR, y =mean, color = FemState ))+
geom_line(linewidth=1, show.legend = NA)+
theme_bw(base_size = 8)+
geom_ribbon(aes(ymax = upr, ymin = lwr, x=OSR, fill=FemState),color=NA, alpha = 0.1)+
ylab(label = "Log-odds ratio of FP vs alternative")+
xlab(label="Log OSR")+ #xlim(0,80)+
scale_color_manual(breaks = c("FPvsSM", "FPvsSD"),
values = c(pal[c(5,9)]),
labels = c("SM", "SD"),
name = "Alternative state") +
scale_fill_manual(breaks = c("FPvsSM", "FPvsSD"),
values = c(pal[c(5,9)]),
labels = c("SM", "SD"),
name = "Alternative state") +
geom_hline(yintercept = 0, colour = "gray50", lty = 2) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
pm4
Hypothesis test
n <- nrow(mm4.1$Sol)
Model4.summ <- data.frame(Log_odds = c("FP_vs_SM", "FP_vs_SD"),
PM = c(mean(mm4.1$Sol[,3]), mean(mm4.2$Sol[,4])),
Lwr = c(HPDinterval(mm4.1$Sol[,3])[1],
HPDinterval(mm4.2$Sol[,4])[1]),
Upr = c(HPDinterval(mm4.1$Sol[,3])[2],
HPDinterval(mm4.2$Sol[,4])[2]),
PMCMC = c(sum(mm4.1$Sol[,3] < 0) / n,
sum(mm4.2$Sol[,4] < 0) / n))
Model4.summ %>%
kbl(booktabs =T, align= c(rep('l',1),rep('r', 4))) %>%
kable_styling(latex_options ="striped")| Log_odds | PM | Lwr | Upr | PMCMC |
|---|---|---|---|---|
| FP_vs_SM | 0.0435336 | -0.8805157 | 0.9239611 | 0.4643333 |
| FP_vs_SD | -0.5332423 | -1.5030749 | 0.3680942 | 0.8750000 |
Is population density higher in the temperate zone?
See introduction in Willink et al. (2024) for the rationale of this hypothesis.
Data summary
First plot raw data
ecodat2$combined <- paste(ecodat2$Lat,ecodat2$Hab, sep = "-")
pdem_eco <- ggplot(data = ecodat2, aes(x = combined, y = log(dens), colour = FemState)) +
geom_jitter(width = 0.2, size = 1.5, alpha = 0.7)+
theme_bw(base_size = 8) +
labs(x = "Latitude-habitat", y = "Log adult density") +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
scale_fill_manual(
breaks = c(0,1,2),
values = pal[c(1,12,7)],
labels = c("SM", "FP", "SD"),
name = "Colour state"
) +
scale_colour_manual(
breaks = c(0,1,2),
values = pal[c(1,12,7)],
labels = c("SM", "FP", "SD"),
name = "Colour state"
) +
scale_x_discrete(labels = c("temperate-closed", "temperate-open", "tropical-closed", "tropical-open")) +
theme(axis.text.x = element_text(angle = 10, vjust = 0.85))
#ggsave(pdem_eco, filename = "~/Dropbox/Doctorado/COEN_Comparative/Writing/Figures/dens_eco1.pdf", device = "pdf", units = "mm", width = 120, height = 80)Model 5: total adults by habitat and sampling effort
In Model 5 we estimated the number of adults for each latitude-habitat combination at the average sampling effort (i.e. by using the mean-centered catching effort as a covariate).
Run the model
inv.phylo<-inverseA(ecotre2,nodes="TIPS",scale=TRUE)
#set prior
prpois1 = list(R=list(V=1,nu = 0.002),
G =list(G1=list(V=diag(1),nu=2,alpha.mu=rep(0,1),alpha.V=diag(1)*1000),
G2=list(V=diag(1),nu=1000,alpha.mu=rep(0,1),alpha.V=diag(1))))
#The model for the interaction term
#starting value
mm5 <- MCMCglmm(Total ~ Hab:Lat + Hab:Lat:TimeC -1,
random = ~Site+Taxon,
ginverse=list(Taxon=inv.phylo$Ainv),
family ="poisson", data = ecodat2,
prior = prpois1,pl=TRUE,
nitt=330000, burnin=30000, thin=100,verbose = F)Plot predictions
#colnames(ecodat2)
newdat <- ecodat2[,c(4,5,6,7,8,11,12,13)]
Pred <- predict.MCMCglmm(mm5, newdata = newdat,interval = "confidence", type = "response")
Pdata <- cbind(newdat, Pred)
Pdata$Eco <- paste(Pdata$Lat, Pdata$Hab, sep = "_")
Pdata$Lat <- recode(Pdata$Lat, "Trp" = "Tropical", "Tmp" = "Temperate")
Pdata <- Pdata[!Pdata$Taxon == "Xanthocnemis_zealandica",]
pal <- wes_palette("Zissou1", 12, type = "continuous")
p5 <- ggplot(data = Pdata, aes(x=TimeC, y =fit, color = Eco ))+
labs(title = "", x="Catching effort", y="Number of individuals") +
geom_line(linewidth=0.5, show.legend = NA)+
# xlim(0,400)+
#ylim(0,80)+
theme_bw(base_size = 8)+
#facet_wrap(~Lat, ncol = 2)+
geom_ribbon(data = Pdata,aes(ymax = upr, ymin = lwr, x=TimeC, fill=Eco),color=NA, alpha = 0.1)+
geom_point(aes(y =Total, fill = Eco), size=1 , alpha = 0.3)+
scale_fill_manual(
breaks = c("Trp_C", "Trp_O", "Tmp_C", "Tmp_O"),
values = pal[c(12, 8, 5, 1)],
labels = c(
"tropical-closed",
"tropical-open",
"temperate-closed",
"t
emperate-open"
),
name = "Latitude-habitat"
) +
scale_colour_manual(
breaks = c("Trp_C", "Trp_O", "Tmp_C", "Tmp_O"),
values = pal[c(12, 8, 5, 1)],
labels = c(
"tropical-closed",
"tropical-open",
"temperate-closed",
"temperate-open"
),
name = "Latitude-habitat"
) +
guides(fill=guide_legend(nrow=2,byrow=TRUE)) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
theme(legend.position = "top")
p5
Hypotheses testing
# Temperate vs tropical in open
LatO_m <- mean(mm5$Sol[,2] - mm5$Sol[,4])
LatO_l <- HPDinterval(mm5$Sol[,2] - mm5$Sol[,4])[1]
LatO_u <- HPDinterval(mm5$Sol[,2] - mm5$Sol[,4])[2]
LatO_p <- sum(mm5$Sol[,2] - mm5$Sol[,4] < 0) / n
# Temperate vs tropical in closed
LatC_m <- mean(mm5$Sol[,1] - mm5$Sol[,3])
LatC_l <- HPDinterval(mm5$Sol[,1] - mm5$Sol[,3])[1]
LatC_u <- HPDinterval(mm5$Sol[,1] - mm5$Sol[,3])[2]
LatC_p <- sum(mm5$Sol[,1] - mm5$Sol[,3] < 0) / n
# Open vs closed in Temperate
HabTmp_m <- mean(mm5$Sol[,2] - mm5$Sol[,1])
HabTmp_l <- HPDinterval(mm5$Sol[,2] - mm5$Sol[,1])[1]
HabTmp_u <- HPDinterval(mm5$Sol[,2] - mm5$Sol[,1])[2]
HabTmp_p <- sum(mm5$Sol[,2] - mm5$Sol[,1] < 0) / n
# Open vs closed in Tropical
HabTrp_m <- mean(mm5$Sol[,4] - mm5$Sol[,3])
HabTrp_l <- HPDinterval(mm5$Sol[,4] - mm5$Sol[,3])[1]
HabTrp_u <- HPDinterval(mm5$Sol[,4] - mm5$Sol[,3])[2]
HabTrp_p <- sum(mm5$Sol[,4] - mm5$Sol[,3] < 0) / n
# combine in dataframe
Dens_tests <- data.frame(Comparsion = c("temperate-open vs tropical-open",
"temperate-closed vs tropical-closed",
"tropical-open vs tropical-closed",
"temperate-open vs temperate-closed"),
Mean = c(LatO_m, LatC_m, HabTr_m, HabTmp_m),
Lower = c(LatO_l, LatC_l, HabTrp_l, HabTmp_l),
Upper = c(LatO_u, LatC_u, HabTrp_u, HabTmp_u),
PMCMC = c(LatO_p, LatC_p, HabTrp_p, HabTmp_p))
Dens_tests %>%
kbl(booktabs =T, align= c(rep('l',1),rep('r', 4))) %>%
kable_styling(latex_options ="striped")| Comparsion | Mean | Lower | Upper | PMCMC |
|---|---|---|---|---|
| temperate-open vs tropical-open | 0.9221925 | 0.3999896 | 1.392118 | 0.0000000 |
| temperate-closed vs tropical-closed | 1.1386435 | 0.3579164 | 1.940382 | 0.0036667 |
| tropical-open vs tropical-closed | 0.0625325 | 0.3574037 | 1.362952 | 0.0003333 |
| temperate-open vs temperate-closed | 0.6593395 | 0.0168744 | 1.309482 | 0.0216667 |
Are ecological effects mediated by demographic effects?
In Model 6 we repeat the analysis of ecological effects on the occurrence of FP (Model 2), but statistically controlling for the demographic effect of adult density. If ecological effects are indeed entirely mediated by density we expect that once we know the density of a population, knowing its habitat and latitudinal range should not change our believe of the probability of FPs. In other words, the large effect of open habitats in temperate regions should disappear once we control for density.
Run the model
# we need a matrix of phylogenetic distances to account for shared evolutionary history
inv.phylo <- inverseA(ecotre2, nodes = "TIPS", scale = TRUE)
# residual covariance matrix 1/J*(I+J), I and J are identity and unit matrices respectively
IJ <- (1/3) * (diag(2) + matrix(1, 2, 2))
# kronecker prior for fixed effect close to being flat for the 2 way:
#(0 vs.1 , 1 vs. 2 and 0 vs. 2) marginal probabilities within each fixed effect
pr6 <- list(
B = list(mu = rep(0, 16), V = diag(16) * (1.7 + pi ^ 2 / 3)),
R = list(V = IJ, fix = 1),
G = list(
G1 = list(
V = IJ,
nu = 1000,
alpha.mu = rep(0, 2),
alpha.V = diag(2)
),
G2 = list(
V = 1,
nu = 1000,
alpha.mu = rep(0, 1),
alpha.V = diag(1)
)
)
)
#Run the model without intercept to estimate the probability of each state in each combination of latitude and habitat
#starting value
mm6.1 <-
MCMCglmm(
relevel(FemState, ref="0") ~ Lat:Hab:trait + Lat:Hab:log(dens):trait -1,
random = ~ us(trait):Taxon + Site,
rcov = ~ us(trait):units,
ginverse = list(Taxon = inv.phylo$Ainv),
family = "categorical",
data = ecodat2,
prior = pr6,
pl = TRUE,
slice = TRUE,
nitt = 1100000,
burnin = 100000,
thin = 500,
verbose = F
)
mm6.2 <-
MCMCglmm(
relevel(FemState, ref="2") ~ Lat:Hab:trait + Lat:Hab:log(dens):trait -1,
random = ~ us(trait):Taxon + Site,
rcov = ~ us(trait):units,
ginverse = list(Taxon = inv.phylo$Ainv),
family = "categorical",
data = ecodat2,
prior = pr6,
pl = TRUE,
slice = TRUE,
nitt = 1100000,
burnin = 100000,
thin = 500,
verbose = F
)Make log-odds dataset
LO_dat <- data.frame(model = c(rep("0", 8), rep("1", 8)),
Lat = rep(rep(c("Tmp", "Trp"), each = 4), 2),
Hab = c(rep(rep(c("C", "O"), each = 1),2), rep(rep(c("C", "O"), each = 1),2)),
comp = rep(rep(c("FP vs SM", "FP vs SD"), each = 2), 8),
mean = c(mean(mm2.1$Sol[,7]),
mean(mm2.1$Sol[,3]),
mean(mm2.2$Sol[,8]),
mean(mm2.2$Sol[,4]),
mean(mm2.1$Sol[,5]),
mean(mm2.1$Sol[,1]),
mean(mm2.2$Sol[,6]),
mean(mm2.2$Sol[,2]),
mean(mm6.1$Sol[,1]),
mean(mm6.1$Sol[,3]),
mean(mm6.2$Sol[,5]),
mean(mm6.2$Sol[,7]),
mean(mm6.1$Sol[,2]),
mean(mm6.1$Sol[,4]),
mean(mm6.2$Sol[,6]),
mean(mm6.2$Sol[,9])),
lwr = c(HPDinterval(mcmc(mm2.1$Sol[,7 ]))[1],
HPDinterval(mcmc(mm2.1$Sol[,3 ]))[1],
HPDinterval(mcmc(mm2.2$Sol[,8 ]))[1],
HPDinterval(mcmc(mm2.2$Sol[,4 ]))[1],
HPDinterval(mcmc(mm2.1$Sol[,5 ]))[1],
HPDinterval(mcmc(mm2.1$Sol[,1 ]))[1],
HPDinterval(mcmc(mm2.2$Sol[,6 ]))[1],
HPDinterval(mcmc(mm2.2$Sol[,2 ]))[1],
HPDinterval(mcmc(mm6.1$Sol[,1 ]))[1],
HPDinterval(mcmc(mm6.1$Sol[,3 ]))[1],
HPDinterval(mcmc(mm6.2$Sol[,5 ]))[1],
HPDinterval(mcmc(mm6.2$Sol[,7 ]))[1],
HPDinterval(mcmc(mm6.1$Sol[,2 ]))[1],
HPDinterval(mcmc(mm6.1$Sol[,4 ]))[1],
HPDinterval(mcmc(mm6.2$Sol[,6 ]))[1],
HPDinterval(mcmc(mm6.2$Sol[,9 ]))[1]),
upr = c(HPDinterval(mcmc(mm2.1$Sol[,7 ]))[2],
HPDinterval(mcmc(mm2.1$Sol[,3 ]))[2],
HPDinterval(mcmc(mm2.2$Sol[,8 ]))[2],
HPDinterval(mcmc(mm2.2$Sol[,4 ]))[2],
HPDinterval(mcmc(mm2.1$Sol[,5 ]))[2],
HPDinterval(mcmc(mm2.1$Sol[,1 ]))[2],
HPDinterval(mcmc(mm2.2$Sol[,6 ]))[2],
HPDinterval(mcmc(mm2.2$Sol[,2 ]))[2],
HPDinterval(mcmc(mm6.1$Sol[,1 ]))[2],
HPDinterval(mcmc(mm6.1$Sol[,3 ]))[2],
HPDinterval(mcmc(mm6.2$Sol[,5 ]))[2],
HPDinterval(mcmc(mm6.2$Sol[,7 ]))[2],
HPDinterval(mcmc(mm6.1$Sol[,2 ]))[2],
HPDinterval(mcmc(mm6.1$Sol[,4 ]))[2],
HPDinterval(mcmc(mm6.2$Sol[,6 ]))[2],
HPDinterval(mcmc(mm6.2$Sol[,9 ]))[2])
)Plot final log odds data
LO_dat$eco <- paste(LO_dat$Lat, LO_dat$Hab, sep = "-")
p6 <- ggplot(data = LO_dat, aes(x = eco, y = mean, colour = model)) +
geom_pointrange(aes(ymin = lwr, ymax = upr), position = position_dodge(width = 0.3), size = 0.4) +
labs(y="Log-odds of FP vs alternative state", x = "Latitude-habitat")+
facet_wrap(~comp) +
geom_hline(yintercept = 0, colour = "grey50", lty=2)+
scale_color_manual(breaks = c("0", "1"),
values = c(pal[c(5,9)]),
labels = c("no", "yes"),
name = "Controlling \nfor density") +
scale_x_discrete(labels = c(
"Temperate\nclosed",
"Temperate\nopen",
"Tropical\nclosed",
"Tropical\nopen"
)) +
theme_bw(base_size = 8) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
theme(axis.text.x = element_text(angle = 15, vjust = 0.85))
p6
#ggsave(p6, filename = "~/Dropbox/Doctorado/COEN_Comparative/Writing/Figures/Log-odds_mediator.pdf", device = "pdf", units = "mm", width =100, height = 70)
Hadfield J.D. 2010. MCMC methods for
multi-response generalized linear mixed models: The
MCMCglmm R package. Journal of Statistical
Software. 33:1–22.
Höhna S., Landis M.J., Heath T.A., Boussau B., Lartillot N., Moore B.R.,
Huelsenbeck J.P., Ronquist F. 2016. RevBayes:
Bayesian phylogenetic inference using graphical models and
an interactive model-specification language. Systematic Biology.
65:726–736.
Landis M., Edwards E.J., Donoghue M.J. 2021. Modeling phylogenetic biome
shifts on a planet with a past. Systematic Biology. 70:86–107.
Willink B., Ware J., Svensson E.I. 2024. Tropical origin, global
diversification and dispersal in the pond damselflies
(Coenagrionoidea) revealed by a new molecular phylogeny.
Systematic Biology. syae004.